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ABSTRACT 


The  non-linear  elasto-pla  Stic  responses  of  a  submerged 
cylindrical  shell  to  an  underwater  shock  wave  have  been 
investigated.  Using  the  EPSA  (Elasto-Plast ic  Shell 

Analysis)  code,  the  gross  responses  of  homogeneous  and 
ring-stiffened  shells  were  evaluated.  The  relevant  parame¬ 
ters  have  been  displayed  and  evaluated  using  PATRAN-G  color 
graphics  system.  An  interface  module  was  developed  between 
EPSA  and  PATRAN-G  .  The  deformations  and  von  Hises  stresses 
throughout  the  shell  have  been  qualitatively  evaluated. 
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I.  IS TR0DDCTI3B 


A.  GENERAL  PRESENTATION  OP  THE  PHENOMENA 

Great  progress  has  been  made  in  the  static  and  dynamic 
analysis  of  complex  structures  through  the  continued  devel¬ 
opment  of  discrete  element  methods  of  structural  analysis. 
Tremendous  improvements  in  computing  power  ha"e  made 
possible  the  study  of  fully  nonlinear  problems. 

The  analysis  of  the  response  of  a  structure  submerged  in 
a  fluid,  is  severely  complicated  by  the  intrusion  of  signif¬ 
icant  fluid-structure  interaction  effects.  Recently,  the 
development  of  a  variety  of  surface  interaction  approxima¬ 
tions  has  provided  the  means  for  a  more  efficient  analysis 
of  the  coupling  between  the  structure  and  the  surrounding 
fluid. 

Computer  codes  for  structural  analysis  are  well- 
developed  so  that  the  fluid-structure  interaction  is,  for 
the  most  part,  handled  by  adding  new  capabilities  to 
existing  structural  analysis  programs. 

It  is  a  well  known  fact  that  the  primary  threat  to  a 
submerged  structure  is  the  shock  wave  that  results  from  an 
underwater  explosion.  However,  the  complexity  of  the  phys¬ 
ical  phenomena  involved,  along  with  the  difficulties  encoun¬ 
tered  in  obtaining  experimental  results  have  been  serious 
drawbacks  for  the  analysis  of  these  types  of  problems.  But 
there  is  a  definite  need  for  investigations  of  large  defor¬ 
mations  and  buckling  problems  in  a  structure  submitted  to  an 
underwater  explosion. 
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B.  OBJECTIVES 


This  study  deals  with  the  nonlinear  response  of  a 
submerged  cylindrical  shell  to  a  shock  wave.  The  existing 
finite  element  code  EPSA  (Elasto- Plastic  Shell  Analysis) 
[Ref.  1]  which  includes  nonlinear  affects  and  a  surface 
interaction  approximation  was  selected  for  the  study.  In 
order  to  alleviate  the  tedious  interpretation  of  results  at 
points  throughout  the  shell,  PATS AN-G,  a  color  graphics 
system,  was  used.  PAIR AN- 3  allows  for  a  global  representa¬ 
tion  of  a  quantity  distribution  over  the  structure  rather 
than  the  discrete  representation  given  by  a  computer  output. 
The  objectives  of  this  study  were  to  merge  the  finite 
element  code  EPSA  with  th9  color  graphics  system  PATRAN-G, 
and  to  conduct  an  analysis  of  the  response  of  various  types 
of  cylindrical  shells  to  a  spherical  shock  wave  generated  by 
an  underwater  explosion. 


II.  SUE  EPS A  COBPPTBB  P50GBAH 


A.  PBESEHT  ATIOI 

EPSA  (Elasto-Plastic  Shell  Analysis)  [Bef.  1]  is  a 
computer  program  developed  by  Weidlinger  Associates  and 
funded  by  DNA/N AVSSA/CNB  for  the  purpose  of  the  analysis  of 
submerged  stiffened  shells  under  dynamic  loadings.  It  incor¬ 
porates  a  number  of  specific  features  which  are  geared  for  a 
more  efficient  analysis. 

In  particular,  EPSA  allows: 

-The  analysis  Of  shells  in  an  acoustic  medium,  subjected  to 
both  low  and  high  frequency  shock  loadings. 

-An  efficient  modelling  of  the  elasto-plastic  behavior 
-The  inclusion  of  large  displacement  effects  to  analyze 
dynamic  buckling  situations  and  post-buckling  behavior. 

-The  modelling  of  stiffeners  and  internal  structures. 

-The  fluid-structure  interaction  effect 

The  following  sections  describe  the  equations  of  motion 
for  a  submerged  structural  shell  in  an  infinite  fluid 
sujected  to  a  pressure  loading  and  investigate  the  modelling 
of  the  surrounding  fluid.  The  last  sections  are  devoted  to 
the  finite  element  discretization  as  well  as  to  specific 
features  concerning  EPSA. 

B.  EQOATIOIS  OF  BOTIOH  FOB  THE  SHELL 

Considering  a  thin  shell  of  thickness  h  ,  volume  V,  area 
B  submerged  in  an  infinite  fluid  (figure  2.1).  The  shell 
stress  resultants  are  defined  from  the  stress  tensor  by  : 
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Figure  2.1  Shell  Stress  State 


Applying  the  principle  of  virtual  work  gives: 


/. 


(s)  {  Se)  dR 


(p}T {6U}  dR 


P  (aHSu}  SR 


(2.1) 


Where 

{u}»(ulru2#w  )T  is  the  displaceaent  vector  at  each  point 


(s>  *  (Nil, N  22 

vector 

*  N  12#H  ii 

22»  H  12  )T 

is 

the 

stress 

resultant 

(e|  *  (en»e22 
vector 

,k  -  2k,  _)T 

2  2  12 

is 

the 

strain 

resultant 

p is  the  mass  per  unit  area  for  the  shell. 

The  symbol  5  will  refer  to  a  virtual  quantity,  and  the  dot 
or  star  denotes  a  differentiation  with  respect  to  time. 


■V 


Tha  first  term  of  equation  (2.1)  rapresents  the  virtual  work 
of  internal  forces  ,  the  second  represents  the  virtual  work 
of  external  forces  (i.e.  pressure,  point  loading,  etc)  ,  and 
the  third  one  represents  the  contribution  of  inertia  forces 
in  the  virtual  work.  Thus,  this  last  term  expresses  the 
effects  of  dynamic  phenomena  on  the  structure. 

C.  FLO  ID  BODE  LING 

In  the  case  of  a  submerged  structure,  the  external 
forces  are  the  pressures  applied  at  the  fluid-structure 
interface. 

As  the  shock  wave  hits  the  stucture  it  gets  reflected, 
thus  inducing  a  pressure  tarm  pr  .  In  addition,  the  motions 
of  the  shell  will  also  generate  a  radiated  wave,  with  a 
pressure  contribution  pra 

Therefore,  the  pressure  at  tha  fluid  structure  interface  is 
the  sum  of  the  incident,  reflected  and  radiated  pressures: 


p  =  Pi  ♦  pr  ♦  Pra  (2.2) 

Where 

p  =  Total  dynamic  pressure . 

P;  =  Pressure  associated  with  incident  free  fiald  pressure 
wave. 

pr  =  Reflected  pressure  due  to  th9  interaction  of  the  inci¬ 
dent  wave  with  the  structure,  the  structure  being  fixed  and 
rigid. 

Pra=  Radiated  pressure  due  to  the  normal  movements  of  the 

surface  of  the  structure  in  the  fluid 

Psa  Pr*  Pra^s  called  the  scattered  pressure. 

The  methods  for  getting  the  scattered  pressure  pg  will 
now  be  investigated.  Assuming  a  spherical  wave  in  an 
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acoustic  medium  with  a  sound  speed  c  ,  the  pressure  is 
determined  by  the  well  known  wave  equation: 


__2 

c*V  P  = 


(2.3) 


and  the  proper  conditions  at  the  boundaries  of  the  fluid 
domain  . 

One  alternative  for  solving  the  previous  problem  would 
be  simply  to  use  a  finite  element  discretization  of  part  of 
the  fluid  domain  ,  imposing  a  radiation  condition  at  the 
boundary  [Ref.  2].  However,  this  would  require  a  very 
significant  fraction  of  the  computing  effort  that  could  not 
be  devoted  to  the  structural  modelling. 

Therefore,  a  more  efficient  way  of  computing  the  scattered 
pressure  is  required. 

The  Doubly  Asymptotic  Approximation  (DAA)  imparts  upon 
the  structural  model  a  surface  loading  composed  of  incident 
and  scattered  waves. 

In  the  high  frequency  limit,  it  can  be  shown  that  the 
scattered  nodal  force  vector  (Fs)  is  related  to  the  wave 
particle  velocities  normal  to  the  structure's  surface  by 

[Ref.  3]  : 


{Fs}  *  [  A  ]C<JS)  (2.4) 

Where  [ A  ]  is  the  diagonal  matrix  of  nodal  areas  (areas  asso¬ 
ciated  with  each  node)  and  (0S)  is  the  vector  of  nodal  scat¬ 
tered  normal  velocities.  Therefore,  in  this  high  frequency 
case  the  shock  wave  is  simply  a  plane  wave  and  equation 
(2.4)  states  that  the  pressure  is  proportional  to  the  fluid 
velocity. 


16 


In  the  low  frequency  Unit  tha  fluid  structure  interac¬ 
tion  is  governed  by  the  relation: 

(Fs)  *  Cav](U*}  <2. 5) 

Share  {0„*  }  =  is  the  nodal  normal  acceleration  vector 

and  [  ]  is  the  added  mass  matrix  computed  in  an  incompres¬ 

sible  fluid. 

Thus,  in  the  low  frequency  case,  the  loading  is  due  to 
tha  motion  of  the  rigid  structure  in  an  incompressible 
fluid,  a  problem  typically  found  in  hydrodynamics. 

When  a  broader  range  of  frequencies  is  considered,  one 
can  combine  the  two  previous  equations  wirh  the  differential 
equation  governing  the  scattered  prassure  [Ref.  3],  giving: 


[A}*  CFS*}  *  {Fsl  *  PcC^*}  (2-6) 

Where  {Fs*}  *  ^  (Fs  } 

Defining  the  vector  cf  nodal  scattered  pressures  (pg)  by: 

CPS]  «  [A]-1  (Fs) 

we  get: 

C«V](PS}  ♦  pc[A](Ps)  «  pc[Hv][0s)  (2.7) 

which  is  the  set  of  equations  that  gives  the  scattered  pres¬ 
sure  at  each  node  of  the  wetted  surface  of  the  shell. 
Equation  (2.7)  gives  exact  results  in  both  the  high  and  low 
frequency  limits.  Thus,  DA A  allows  the  modelling  of  the 
fluid-structure  interaction  via  a  coupled  set  of  differen¬ 
tial  equations  at  the  wetted  nodes  of  the  structure  instead 
of  modelling  the  whole  fluid  with  a  finite  element  grid. 
The  use  of  DAA  will  free  some  memory  space  in  the  computer 


for  a  better  modelling  of  structural  behavior  while  at  the 
same  time  giving  soae  reasonably  accurate  results  for  the 
loading  of  the  structure. 

Therefore,  the  equations  for  the  study  of  underwater 
shocks  will  consist  of  a  coupled  set  of  structural  equations 
that  come  froa  the  Principle  of  Virtual  Work  and  of  fluid 
equations  at  the  wetted  nodes  that  roae  from  DA  A  equations. 


D.  FIHITE  ELEMENT  PfiOCEDOHE 

1 .  Discretization 


The  principle  of  virtual  work  is  rewritten  for  the 
structure  introduced  in  section  B  [Bef.  4] 


Cde}  dS  - 


(5u}dR  ♦ 


u}dR  =  0 


(2.8) 


The  surface  of  the  region  is  covered  by  a  quadrilat¬ 
eral  mesh  cf  N  elements,'  each  having  an  area  Ai  .  The 
previous  integral  can  than  be  replaced  by  a  summation  of 
integrals  over  A ^  .  : 


{  Se}  d  R  = 


(5e)i  dR 


(2.9) 


The  integrations  over  A .  ar e  then  performed  by  dividing  the 

^  k 

element  into  four  regions  AA  (figure  2.2)  He  have  then: 


(s (5e).  dR 


t 

k=: 


(s}iT  (Se*^ 


k 


(2. 10) 


and  therefore  equation  (2.9)  becomes  : 


{ 6  e}  d  R 


-  t.  (3li  tS9kii  f 

i* 1  k” 1 


Using  the  same  procedure,  it  can  also  be  derived  : 


(2.11) 
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f  T  ^  A'  T  „  k  T  „ 

/  { P }  (5u)  dR=  2^2^  tPi  €5  A±  =  {P}  {6q} 

■'R  i  =  l  k=  1 


f  ..  N  4 

/p{u}T(6u}  dR=  £  £  CuJi  {du}i  4=  C M 3  Cq}  c«q> 

"k  i  =  l  k=l 


(2.  12) 


(2.13) 


Where  [  M  ]  is  the  aass  matrix  ,  {?}  is  the  vector  of  exter¬ 

nally  applied  forces,  and  {q}  is  the  vector  of  nodal  normal 
displacements  for  the  structure. 


Figure  2.2  Srid  Points  in  EPS&. 

By  definition,  finite  element  discretization  can 
express  the  displacement  {u}  at  any  point  of  an  element  as  a 
function  of  the  displacements  at  the  corner  points  of  the 
element,  defined  by  {q}  . 

(U)  *  [H  ]  {q}  (2.  14) 


Where  [H]  is  a  6x12  matrix  of  interpolation  functions. 

Combining  the  derivatives  of  (uj  will  give  the  strain  vector 
(e) .  In  matrix  form,  equation  (2.14)  gives  : 


{9} 


[B]  {<?} 


(2.  15) 


Where  [B]  is  a  6^12  matrix  function  of  the  geometry  of  the 
element  as  well  as  a  function  of  the  previous  deformations 
in  the  case  of  large  displacements. 

Using  the  previous  result,  equation  (2.9)  is  rewritten  in 
the  following  way  : 


f{s}T(6e}  da*  £  £  {s}J‘[Bk]i{5q}i  {P}T  {6q} 

i= lk= 1 

Where  (F)  is  the  internal  force  vector. 


(2.16) 


Combining  the  previous  equations  in  equation  (2.8) 
principle  of  virtual  work  becomes: 


,  the 


(2.  17) 


Therefore,  the  principle  of  virtual  work  has  been  tranformed 
into  a  system  of  ordinary  differential  equations  which  are 
much  more  amenable  to  numerical  treatment. 


In  EPS  & ,  each  arbitrarily  shaped  quadrilateral 
element  is  defined  by  four  corner  nodes,  each  having  three 
translational  and  no  rotational  degree  of  freedom.  In  order 
to  represent  the  behavior  in  bending  eight  nodes  not  contig¬ 
uous  with  the  element  are  also  used  (figure  2.3).  Every 
element  accesses  twelve  nodes  and  has  twenty  degrees  of 
freedom:  three  t ranslati cnal  degrees  at  the  corner  nodes 
and  one  corresponding  to  the  displacement  normal  to  the 
surface  for  each  of  the  eight  exterior  nodes.  The  bending 
behavior  (second  derivative  term)  is  expressed  in  terms  of 
the  nodal  displacements  via  a  finite  difference  technique 
CRef.  «]. 
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Then  the  nodal  dis placement  vector  of  element  i  is 


simply: 


{g}  =  <uira2  ,u;j.a4,v1,. .  .v4,w1 .w4  ,v5  .w12  T  (2.18) 

Conventional  finite  element  codes  utilize  three 
translational  and  two  rotational  degrees  of  freedom  at  each 
node  (each  element  has  5  4  =  20  d.o.f.  as  in  EPSA  ).  However 
the  masses  associated  with  rotational  degrees  of  freedom  are 
very  small,  leading  to  numerical  instability.  The  use  of 
the  aforementioned  formulation  alleviates  this  problem  since 
only  translations  are  considered  and,  in  addition,  the 
number  of  unknowns  is  reduced,  leading  to  simpler  and  more 
efficient  computations. 

It  must  be  pointed  out  that  in  order  to  model  the 
bending  behavior  at  the  edges  of  the  shell,  a  set  of  artifi¬ 
cial  nodes  has  to  be  created.  The  finite  element  grid  will 
then  consist  of  the  nodes  defined  in  the  input  file  plus  the 
artificial  nodes  along  the  boundary  of  the  sheet 
(figure  2.3). 

2 .  Strain  Displacement  Relation 

The  principle  of  virtual  work  deals  only  with 
strains.  Since  the  finite  element  approximation  gives  the 
displacement  at  each  point,  the  equations  relating  the 
strains'  to  the  displacements  are  needed.  In  this  study,  the 
Donnell- Vlasov  nonlinear  kinematic  equations  of  shell  theory 
are  employed,  and  the  strain-displacement  relations  are 
described  in  table  I. 

Using  equation  (2.15)  in  the  equations  detailed  in 
the  previous  section  will  give  the  finite  element  approxima¬ 
tion  of  strains  in  terms  of  nodal  displacements. 
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T  ABLE  I 

Donne 11-flasov  Shell  Equations 
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l  '£  are  local  coordinates 
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h1  ,h2  are  scaling  coefficients. 
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3.  ihfiii  £an§titu£iTe  Relations 

The  shell  constitutive  relations  relate  the  stress 
resultant  rate  vector  to  the  shell  strain  rate  vector.  In 
natrix  terns: 


An 


Figure  2.3  Nodal  Points  Drganxzation. 


The  stress-strain  relation  used  in  EPSA  differs  from 
the  classical  Elast o-Plastic  theory  in  that  the  formulation 
involves  shell  stress  resultants  rather  than  stresses  at 
points  throughout  the  thickness  of  the  shell.  In  other 
words,  EPSA  uses  relation  2.19  integrated  over  the  thickness 
of  the  shell.  Thus,  there  is  no  need  to  compute  the  stresses 
throughout  the  shell,  which  results  in  a  significant  savings 
in  storage  space  and  processing  time.  However,  the  stress 
resultants  N^and  cf  th9  shell  theory  are  not  sufficient 
to  describe  the  state  of  stress,  and  certain  higher  order 
moments  must  be  combined  with  the  stress  resultants  to  form 
the  dynamic  variables  of  the  problem.  The  coefficients  for 
the  integrated  constitutive  equation  have  been  determined 
using  experimental  results  [Bef.  5). 
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These  constitutive  equations  consist  of  a  yield  condition,  a 
strain  hardening  law  and  a  flow  rule: 

-The  yield  condition  indicates  whether  part  of  the  shell  has 
started  to  yield.  (figure  2.4) 

-The  strain  hardening  law  gives  the  evolution  of  stresses  in 
the  shell  after  plasticity  is  reached. 

-The  flow  rule  gives  the  plastic  strain  rate  in  the  shell 
after  plasticity. 


Figure  2.4  Yield  Situation  in  the  Shell. 


4 -  Solution  Procedure 

EPS  A  uses  an  explicit  central  difference  scheme  to 
solve  the  virtual  work  and  fluid  loading  equations  detailed 
in  section  B.  As  indicated  in  appendix  B,  the  explicit  time 
integration  procedure  requires  a  small  time  step  and  is  only 
conditionally  stable.  However  when  stable,  it  always 
converges  to  the  exact  solution,  as  opposed  to  implicit 
schemes  that  are  unconditionally  stable  but  may  lead  to 
unrealistic  rasults.  Furthermore, in  problems  involving  the 
treatment  of  shocks,  accuracy  requirements  preclude  the  use 
of  large  time  steps.  The  central  difference  scheme  seems, 
therefore,  particularly  optimal. 
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Assuming  we  know  the  solution  at  time  z  ,the  central  differ¬ 
ence  scheme  applied  to  equation  (2.17)  gives  the  solution  at 
time  t+At  : 


{7}^=  (V}*  ♦  At  C  CPI  i  -  E  {?} k  )  (2.20) 

1  1 h  k=1 

Where  M  i  is  the  mass  of  node  i  ,  £P)i  are  the  externally 
appied  forces  and  {F)k  are  summed  over  all  the  elements  k 
framing  node  i. 


The  formulation  of  the  equations  is  in  the  initial 
configuration  and  all  equations  are  solved  in  the  initial 
geometry  in  accordance  with  the  total  Lagrangian  formulation 

[Bef.  6]. 


E.  EPSA  CAPABILITIES 

The  structure  to  be  analyzed  is  conceptually  divided 
into  constitutive  parts  called  "sheets."  Each  sheet  is  a 
curved  section  of  a  shell  with  an  arbitrary  number  of  nodes 
and  elements  (figure  2.5)  Its  shape  is  limited  to  a  surface 
that  can  be  described  by  a  smooth  continuous  function 

without  discontinuities  in  its  slope.  The  elements  within 

the  sheet  are  limited  to  a  rectangular  organisation  (figure 
2.5)  . 

Thus,  a  cylinder  with  end  caps  would  consist  of  three 

sheets:  a  cylindrical  sheet  and  a  sheet  for  each  end 

(figure  2.5).  Three  sheets  are  required  because  of  the 
slope  discontinuity  at  the  edge  between  the  cylinder  and  the 
end  caps. 

Dividing  the  structure  into  shaets  isolates  the  diffi¬ 
culties  associated  with  the  boundaries  into  a  set  of  artifi¬ 
cial  nodes  along  the  perimeter  of  the  sheet.  when  several 
sheets  are  required  ,EPSA  takes  care  of  the  compatibilities 
of  displacements,  rotations  and  moments  along  the  edges. 
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Figure  2.5  Sheet  Organization 


Thus,  any  arbitrarily- shaped  structure  can  be  analyzed 
using  EPSA  by  dividing  it  into  a  number  of  sheets. 

Two  options  for  choosing  elements  are  available  in  EPSA. 
The  first  option  exists  to  employ  a  generalized  quadrilat¬ 
eral  element.  The  second  op-cion  exists  to  employ  a  rectan¬ 
gular  element  and  uses  a  specialized  formulation  for  this 
type  of  element. 

The  coordinates  which  can  be  either  cartesian  or  cylin¬ 
drical  always  lie  within  the  plana  of  the  sheet.  The  z 
direction  lies  in  a  positively  outward  direction  normal  to 
the  sheet.  Each  sheet  contains  its  own  local  coordinate 
system,  there  is  no  global  coordinate  system  when  multi 
sheets  are  merged  (figure  2.6). 

Prior  to  generating  a  finite  element  mesh  for  a  sheet 
one  must  establish  the  side  numbers  of  the  sheet.  The  side 
numbering  scheme  is  arbitrary  as  to  the  choice  of  sheet 
number  one.  However  the  specification  of  sides  1  to  «  must 
proceed  in  a  counterclockwise  direction  when  the  sheet  s 
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Sheet  Subdivision 
Coordinate  SvitgB 


ORTHOGONAL  ALONG 
FRIKCIfAL  RAIII 

or  cvxvatiiii 


Figure  2.6  Coordinate  System. 

viewed  from  the  positive  z  direction.  At  the  intersection 
of  side  1  and  side  4,  element  1  is  first  defined.  Then  the 
rode/element  number  is  incremented  by  one  until  it  reaches 
side  2.  Then  it  returns  to  side  4  and  continues  the  count 
for  the  next  row  of  nodes/elements. 

Thanks  to  the  exclusive  use  of  quadrilateral  elements 
and  to  the  specific  numbering  procedure,  the  table  of 
connectivity  is  implicitly  defined  when  generating  the  nodal 
points,  thus  simplifying  considerably  the  input  require¬ 
ments. 

The  inclusion  of  structures  internal  to  the  cylindrical 
sheet  is  enacted  in  EPS  A  through  the  use  of  internal  sheets. 
Structures  internal  to  a  cylinder  are  therefore  modelled  as 
individual  sheets  having  the  same  capabilities  as  any 
general  EPSA  sheet. 

Two  types  of  internal  structures  are  available: 

-Sheets  (beams,  plates  or  shells)  whose  connection  to  the 
cylindrical  shell  lies  parallel  to  the  axis  of  the  cylinder. 
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-Sheets  (beams,  plates,  shells)  whose  connection  lies  in  the 
circumferential  direction  of  the  cylinder. 

Furthermore, in  order  to  modal  internal  equipment 
(machinery,  etc  )  which  does  not  deform  but  contributes  to 
the  nertia  of  the  system,  concentrated  masses  can  be  input 
at  u-dal  points. 

The  user  must  be  awar*  that  the  previously  discussed  DAA 
is  only  implemented  for  cylindrical  stuctures.  Prior  to  the 
finite  element  analysis  the  user  must  compute  the  added  mass 
(virtual  mass)  matrix  defined  in  equation  (2.4).  This  is 
done  by  using  the  ACCESION  program,  which  creates  a  virtual 
mass  file  that  EPSA  reads  when  computing  the  flu  id- structure 
interaction. 

Finally,  EPSA  in  its  latest  version  takes  local  cavita¬ 
tion  into  account.  When  the  total  pressure  computed  by  EPSA 
is  found  to  be  negative,  it  is  simply  set  to  zero  since 
water  cannot  withstand  any  tension. 

F.  USING  EPSA 

The  input  file  for  EPSA  can  be  generated  either 
directly,  or  via  the  interactive  program  PBEPSA  that  prompts 
the  user  to  give  the  input  data  .  Nhen  the  input  file  is 
created,  all  the  data  are  in  free  format. 

The  nodal  points  can  be  generated  semi-automat ically 
(see  the  user's  manual),  and  this  option  is  very  helpful  and 
time  saving  when  generating  big  models.  EPSA  can  be  run 
interactively  for  small  models  or  on  batch  for  bigger  jobs. 
For  instance,  a  cylindrical  shell  with  1440  elements  and 
1517  nodal  points  takes  0.0  129  sec.  of  VAX  CPU  per  time  step 
per  element.  The  whole  model  would  take  about  half  an  hour 
for  100  steps. 

EPSA  creates  an  output  file  in  which  all  input  data  is 
echoed,  and  outputs  the  pressures,  stresses,  strains,  veloc- 


28 


ities,  displacements  at  the  nodal  points  requested  by  the 
user  in  the  input  deck. 

The  value  of  the  time  increment  At  should  be  selected  so 
that  : 

it  <  V2  6nin(  P/E)15  (2.21) 

Whered  .is  the  smallest  distance  between  the  midpoints  of 
opposite  sides  of  an  element,  for  all  elements  of  the  sheet 
(figure  2.5),  and  1/2  is  a  safety  factor.  In  other  words, 
the  time  step  increment  has  to  be  less  than  the  time  taken 
by  a  wave  to  propagate  from  an  element  to  another.  In  equa¬ 
tion  2.21,  (S/p)15  is  simply  the  wave  speed  in  the  material. 

The  virtual  mass  array  (VMA)  is  created  on  unit  20, 
therefore  one  should  not  use  this  unit  for  any  other  purpose 
than  READ  operations. 

Because  of  the  simple  organization  of  its  input  file, 
EPSA  has  been  found  fairly  easy  to  use.  The  user  can  perform 
major  changes  in  the  model  very  quickly  and  efficiently.  The 
ACESION  module  has  been  found  to  work  well  except  for  cylin¬ 
ders  of  large  dimensions  (L=1400",  D=240")  for  which  the 
size  of  the  virtual  mass  array  grows  unexpectedaly  from  a 
reasonable  4  blocks  to  190  blocks  of  VAX  memory  size. 

However,  EPSA  has  been  designed  for  some  specific  types 
of  f luid-stucture  interaction  analysis  and  its  limitations 
should  be  pointed  out: 

-  Only  beam  type  stiffeners  can  be  considered 

The  fluid  structure  interaction  is  only  enacted  for  a 
circular  cylindrical  sheet,  which  takes  away  much  of  the 
flexibility  the  program. 


III.  BPSA/  PATRAN-G  INTERFACE 


A.  INTRODOCTIOH  TO  COLOH  GRAPHICS  SISTBHS 

In  dealing  with  the  response  of  a  structure  to  a 
loading,  quantities  such  as  stresses,  strains,  velocities 
and  displacements  must  be  analyzed  at  a  number  of  nodal 
points,  which  makes  the  interpretation  of  computer  outputs 
very  tedious.  Color  graphics  systems  offar  an  effective 
solution  to  this  problem  by  providing  a  global  representa¬ 
tion  of  a  quantity  distribution  over  a  stucture  rather  than 
a  discrete  representation  given  by  a  computer  output.  A 
color  graphics  system  consists  of  an  interface  between  the 
computer,  the  user  and  the  color  terminal.  A  typical  system 
would  be  a  software  package  that  allows  the  user  to  create 
models  on  the  screen  as  well  as  to  display  any  data  in  terms 
of  color  intensity.  It  is  known  that  a  picture  is  worth 
several  hundreds  of  words.  Therefore,  merging  a  finite 
element  program  with  a  color  graphics  system  would  be  a 
major  improvement  in  engineering  analysis. 

PATRAN-G  [Hef.  7]  is  a  color  graphics  system  specifi¬ 
cally  designed  for  finite  elements,  it  permits  the  engineer 
and  the  computer  to  work  together  towards  the  creation  of  a 
modal.  The  designer  creates  an  image  on  the  screen  and 
PATRAN  automatically  translates  tha  physical  model  into  a 
standard  finite  element  input  deck.  Another  important 
feature  of  PATRAN-G  is  its  ability  to  assist  the  user  in  the 
interpretation  of  results  through  its  post-processing  facil¬ 
ities  which  include:  the  deformed  geometry  with  magnified 
deformations, the  color  coding  of  elements  based  upon  any 
response  parameters  such  as  displacements,  stresses  and 
strains.  In  particular,  the  contour  levels  of  the 
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aforementioned  quantities  can  be  superimposed  on  a 
3-dimensional  image  of  the  model,  thus  allowing  for  a  global 
analysis  of  a  complex  structure. 

B.  HEBGIBG  EPS A  AMD  PATRAN-G 

As  described  in  chapter  I  ,  the  structures  under  study 
have  fairly  short  EPSA  input  files.  Furthermore,  in  this 
study  dealing  with  fluid-structure  interactions  on  a 
circular  cylinder,  only  structures  that  consist  of  one  sheet 
are  considered.  Therefore,  because  of  the  simplicity  of 
both  the  input  file  and  the  stucture  under  study,  there  was 
no  need  to  design  a  module  converting  a  PATRAN-G  model  that 
is  created  on  the  screen  into  an  EPSA  input  file. 

The  remaining  tasks  were  the  following: 

-Display  the  original  finite  element  model  defined  in  the 
EPSA  input  file  on  the  screen  (original  geometry)  . 

-Display  the  nodal  points  and  element  results  that  are 
computed  from  EPSA  on  the  screen  (postprocessing)  . 

1 •  Original  Geometry 

The  input  of  a  finite  element  model  into  PATRAN-G 
can  be  done  by  creating  a  "neutral"  file,  as  described  in 
PATRAN-G  teninology.  The  neutral  file*  is  intended  to 
provide  a  simple  link  between  PATRAN-G  and  the  outside 
world.  It  is  written  entirely  in  80  character  card  images 
and  all  the  data  is  organized  in  small  "packets"  of  two  or 
more  card  images.  Each  packet  contains  the  data  for  a 
fundamental  unit  of  the  model  such  as  node  coordinates  or 
elements  definition.  Since  our  only  purpose  was  to  display 
the  original  geometry  of  the  structure,  a  limited  number  of 


*  Additional  information  about  neutral  files  can  be  found 
in  the  PATRAN-G  user's  manual  [Ref.  7] 


data  packets  (4)  was  to  be  created: 

-File  title  (packet  25)  :  two  cards  containing  the  user 

title. 

-Summary  data  (packet  26)  :  two  cards  containing  the  number 
of  nodes,  elements  and  the  date  and  time  at  which  the 
neutral  file  was  created. 

-Node  data  (packet  1):  contains  all  information  concerning 

nodes:  node  number  and  coordinates  in  a  global  coordinate 

frame. 

-Element  data  (packet  2)  :  contains  the  connectivity  table 

for  the  finite  element  model. 

-End  of  file  (packet  99)  :  end-of-file  cards. 

We  have  seen  in  the  presentation  of  EPSA  in  chapter  I  that 
the  nodes  are  defined  in  a  local,  sheet-attached  coordinate 
system,  that  artificial  nodes  are  created  along  the  edges  of 
the  sheet  to  aodel  the  bending  behavior,  and  that  nc  connec¬ 
tivity  table  wah  input.  Therefore,  the  translator  module 
that  would  be  created  had  to: 

-skip  the  set  of  artificial  nodes  and  properly  renumber  the 
grid 

-perform  a  change  of  coordinates  for  all  local  data 
-generate  the  connectivity  table. 

It  was  decided  to  employ  a  modular  design  in  which 
each  routine  would  perform  a  specific  task.  A  modular  design 
would  allow  further  changes  to  be  made  quickly  and  effi¬ 
ciently  by  modifying  only  the  relevant  routines.  The  imple¬ 
mentation  in  EPSA  was  made  using  a  series  of  "calls",  thus 
minimizing  the  risk  of  interference  with  the  finite  element 
computations. 

The  translator  module  craated  was  made  of  four 
routines: 


-PBELIM  :  computes  the  number  of  nodal  points,  the  number 
of  elements  and  displays  the  first  two  data  packets  (25,26) 

-SHEETF  :  scans  through  the  nodes,  skips  artificial 
nodes,  renumbers  the  grid,  performs  the  required  changes  of 
coordinates  and  displays  the  node  data  packet  (3  1) 

-SHCONN  :  scans  through  the  nodes,  connecting  each  node 
to  the  elements  it  belongs  to  and  displays  the  element  data 
packet  (02)  on  the  neutral  file. 

-ENDNEU  :  displays  the  end-of-file  packet  (99) 

2.  Using  the  Tr anslato r  Module 

The  translator  calls  were  implemented  in  the  routine 
REPORT  of  EPSA.  Any  run  of  EPSA  creates  a  neutral  file  on 
unit  19.  The  neutral  file  name  is  therefore  FOR019.DAT  if  no 
"ASSIGN"  statement  has  been  issued  prior  to  the  computer 
run.  The  finite  element  model  might  than  be  displayed  on  the 
graphic  terminal  (Ramtek,  Tektronix)  via  the  neutral  input 
mode  of  EPSA  (see  [Hef.  7)  for  more  details).  An  example  of 
finite  element  model  output  on  Tektronix  4014  is  given  on 
figure  3.1  . 

3 .  Implementation  of  P ost processing  Capabilities 

Postprocessing  capabilities  exist  to  assist  the  user 
in  the  interpretation  of  computer  results.  It  is  simply  a 
process  of  generating  displays  and  reports  based  upon  a 
combination  of  geometry  and  the  results  of  an  analysis. 

The  results  of  analyses  are  transmitted  to  PATRAN-G 
for  postprocessing  via  "neutral  results  files"  as  decribed 
in  PATRAN-G  terminology.  Unlike  the  model  neutral  file 
described  in  the  previous  section,  results  files  are  in 
binary  rather  than  in  card  image  form. 
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Figure  3.1  Example  of  Pinite  Element  Model  Display 


One  can  distinguish  between  two  major  kinds  of  post¬ 
processing  displays:  deformed  geometries  and  element 
quantities. 

a.  Deformed  Geometry 

A  displacement  results  data  file  required  by 
PATHAN-G  had  to  be  created.  Again,  the  module  created  had  to 
skip  the  artificial  nodes,  perforin  changes  of  coordinates  as 
well  as  to  write  the  nodal  diplacements  in  the  neutral 
results  file.  The  displacement  results  data  file  was 
created  in  module  NEODISP.  Its  organization  is  given  on 
table  II. 

A  small  module  PLOTDISP  that  decides  at  which 
time  steps  the  results  should  be  output  was  created.  The 
"call"  to  PLOTDISP  was  implemented  in  module  COMPOTE  of 
EPSA.  The  overall  structure  of  the  translator  module  is 
presented  on  table  IV  at  the  end  of  the  chapter. 


TABLE  II 

Organisation  of  Displaceman t  Results  File 


!  Column 


Content 


X  displacement  in  global  coordinates 
I  displacement  in  global  coordinates 


Displacement  normal  to  the  shell  without 
rigid  body  mode. 

Velocity  normal  to  the  shell 


A  sample  of  deformed  geometry  output  on  Tektronix  4014  is 
given  in  figure  3.2. 

b.  Element  Quantities  Postprocessing 

Any  element  related  quantity  can  be  the  target 
of  postprocessing.  In  general  these  types  of  quantities  are 
the  results  of  finite  element  computations  supplied  to 
PATRAN-G  through  the  neutral  element  results  file.  The 
neutral  element  results  file  is  different  from  the  neutral 
displacements  results  file  detailed  in  the  previous  section, 
however  it  shares  a  similar  format  ^Bef.  7J.  Element  quan¬ 
tities  such  as  stresses  in  local  and  global  coordinates  and 
von  Hises  stresses  are  computed  in  module  NEtJSTRE  whose 
overall  structure  is  similar  to  NECJDISP  described  earlier. 


Figure  3.2  Deforn9d  Geometry  Output 


Tha  organisation  of  the  neutral  alaaent  results  file  is 
given  on  table  III. 

As  described  in  chapter  I,  EPSA  does  not  compute 
tha  stresses  through  the  thickness  of  the  shell.  Instead  one 
uses  the  integrated  quantities  of  the  shell  theory  [Ref.  8] 


,  W2 

h/2 

r 

/CTij  « 

and 

tdt 

■A/2 

One  can 

expect 

the  stresses  on  the  shell  to  be 

maximum  at  the  extreme  fibers.  NEUSIHE  computes  the  von 
Mises  stresses  at  the  top  and  bottom  fibers  and  writes  the 
maximum  value  in  the  neutral  file.  At  the  surface  of  the 
shell  no  shear  stresses  are  involved.  Assuming  a  linear 
distribution  of  bending  stresses  and  a  uniform  distibution 
of  membrane  stresses,  one  can  easily  derive  : 


°ij  *  Ni-j  /h  ♦  Mij  6/h2  (3.1) 

The  first  term  of  the  previous  equation  is  the  membrane 
force  contribution,  the  second  is  the  bending  moment  contri¬ 
bution.  The  von  (lises  stresses  are  then  computed  using  the 
well-known  relation: 

-t/K  -°1«!  *°2  I  I3-2' 

In  a  similar  way  to  the  displacements  results 
file,  a  module  PLOTSTRE  that  derides  whether  or  not  to 
output  the  element  results  was  created  and  called  from 
COMPOTE.  The  overall  structure  of  tie  translator  is  given  on 
table  IV. 

1 - 

|  T&BLX  III 

j  Organisation  of  the  leotral  Element  Results  Pile 


Column 

CO££2J)t 

Descri  ption 

22 

st re, 1 

Element  local  stress. 

direction  i 

23 

st re, 2 

fl  If  ff 

"  3 

25 

stre, u 

Element  global  stress 

,  direction  x 

26 

stre,  5 

i»  it  ii 

"  y 

27 

stre, 6 

Element  global  shear. 

direction  xy 

31 

von 

von  Hises  stress 

c,  Displaying  the  Results  Quantities  on  the  Screen. 

The  EPS  A  input  deck  has  been  modified  so  as  to 
create  the  results  files  (displacements,  elements)  required 
by  PATHAN-G.  At  the  end  of  the  second  input  card  the  user 
specifies  the  number  of  displacement  results  and  the  number 
of  element  results  files  to  be  created  (at  least  one). 
Obviously,  those  two  inputs  are  also  in  free  format.  The 
neutral  results  files  (elements,  displacements)  will  be 
generated  at  equal  time  intervals  as  the  computation 
proceeds.  The  results  corresponding  to  the  last  time  step 
are  always  output. 

The  element  results  file  is  created  on  unit  16 
and  the  displacement  results  file  on  unit  18,  thus  corre¬ 
sponding  to  files  FOS016.DAT  and  FDR 0  18 .  DAT  respectively.  A 
new  version  of  those  files  is  created  each  time  an  output  is 
requested. 

If  five  (5)  neutral  element  results  files  are  requested  on 
the  input  card  of  a  run  of  20  steps,  five  files  FORO 16. DAT; 1 
to  FOR0 16. DAT; 5  will  be  created,  corresponding  to  time  steps 
u  to  20  respectively. 

For  the  displaying  of  el  ament  and  nodal  points 
results,  the  user  will  refer  to  [Bef.  7]  section  20.  The 
title  of  the  run  (first  card  of  EPSA  input  deck)  will  be 
displayed  on  the  screen  along  with  the  time  at  which  the 
results  have  been  requested. 
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w^csTjaton 


TABLE  IT 

Structure  of  the  Translator 


CALL  PLOTDIS? 


no  output 
requested 


CALL  PLOTS I8E 

no  output 
requested 


computation 

proceeds 


V 


[call  prelimI 
R I  I 

E  I  CALL  SHEET? | 

p|  ! 

OICALL  SHCONNj 
R(  I 

T  CALL  ENDNEUj 


•output  reguested- 


h*— output  requested 


CALL  NEUSTRE 


sea  ns  throug 
the  elements 
change  coord 
ina  tes- 
esmpute  Ton- 
Mises  stress 
-display 
ele  aent 
res  ults 


hi 

-I 

-I 

I 

I 

I 

I 


T 


RETURN 


-CALL  NEUDIS? 


-t- 


scans  throuah  | 
the  grid-  '  I 
changes  coor-  | 
dinates-find  | 
node  with  mar 
defor  mation- 
display  first 
card- 


f 


RETURN 


scans  through  | 
the  grid-  I 

display  nodal  j 
results  I 


I  I 
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IV.  DESCRIPTION  OP  £3E  ORDER RATER  EXPLOSION 


A.  PRESENTATIOH 

An  explosion  is  a  chemical  or  nuclear  reaction  in  a 
substance  (water)  that  converts  an  original  material  into 
other  products  plus  a  significant  amount  of  energy.  The 
process  of  the  explosion  produces  vary  high  temperatures  and 
pressures  and  occurs  with  extreme  rapidity.  As  the  result  of 
the  explosion,  the  initial  mass  of  explosive  becomes  a  very 
hot  mass  of  gas  at  tremendous  pressures;  it  is  then  obvious 
that  these  conditions  will  affect  the  surrounding  medium. 

The  fact  that  the  water  is  compressible  leads  to  the 
conclusion  that  the  pressure  applied  at  some  location  in  the 
liquid  will  propagate  through  it  as  a  wave  disturbance 
[Ref.  9].  It  must  be  pointed  out  that  the  pressures 
involved  in  an  underwater  explosion  are  usually  so  large 
that  the  wave  velocity  cannot  be  assumed  independent  of 
pressure.  Thus,  the  small  amplitude  wave  theory  detailed  in 
[Ref.  9]  does  not  apply  and  this  will  be  the  source  of  many 
complications  in  describing  the  behavior  of  the  shock  wave. 

The  first  cause  cf  disturbance  to  the  water  in  an  under¬ 
water  explosion  is  the  occurrence  of  the  pressure  step  at 
the  water  boundary.  Immediately  upon  its  arrival,  the  pres¬ 
sure  begins  to  be  relieved  by  an  intense  pressure  wave  and 
outward  motion  of  the  water.  For  explosives  like  TNT  or  for 
nuclear  explosions,  the  pressure  rise  can  be  considered  as  a 
discontinuous  step,  and  is  then  followed  by  a  roughly  expo¬ 
nential  decay.  The  duration  of  the  phenomenon  is  of  the 
order  of  a  few  milliseconds  (more  for  nuclear  explosions) . 
Once  initiated,  the  pressure  disturbance  is  propagated  radi¬ 
ally  outward  as  a  compression  wave,  also  called  a  shock  wave 
because  of  the  steep  pressure  step  at  its  front. 


40 


Close  to  the  explosion,  the  velocity  of  the  wave  is  several 
times  the  "acoustic”  speed  of  the  snail  amplitude  theory 
(c=5000  ft/sec)  ;  but  approaches  this  limit  very  rapidly  as 
the  wave  advances  outward. 

The  pressure  level  in  the  wave  falls  more  rapidly  with  the 
radial  distance  than  what  is  predicted  with  the  small  ampli¬ 
tude  theory,  but  approaches  this  behavior  at  large 
distances. 

B.  BOBBIE  EFPECT 

is  a  result  of  the  explosion,  the  initial  mass  of  explo¬ 
sives  becomes  a  hot  mass  of  gas  at  tremendous  pressures. 
After  the  principal  part  of  the  shock  wave  has  been  emitted, 
the  gas  pressure  is  considerably  decreased  but  is  still  much 
higher  than  the  equilibrium  pressure.  The  water  in  the  imme¬ 
diate  region  of  the  sphere  or  "bubble"  of  gas  has  a  large 
outward  velocity  and  the  diameter  of  the  bubble  increases 
rapidly.  The  expansion  continues  and  the  internal  gas  pres¬ 
sure  decreases  gradually,  but  the  motion  persists  because  of 
the  inertia  of  the  outward  flowing  water.  when  the  gas 
pressure  falls  below  the  equilibrium  value,  the  pressure 
defect  brings  the  outward  flow  to  a  stop  and  the  boundary  of 
the  bubble  begins  to  contract  at  an  increasing  rate.  The 
inward  motion  continues  until  the  compressibility  of  the  gas 
reverses  the  motion.  Thus,  the  inertia  of  the  water 
together  with  the  elastic  properties  of  the  gas  provide  the 
necessary  conditions  for  an  oscillatory  system.  The  oscilla¬ 
tions  of  the  gas  sphere  may  persist  a  number  of  cycles,  ten 
or  more  oscillations  having  been  detected  in  favorable 


C.  SURFACE  EPFBCTS 


In  the  case  of  explosions  occurring  at  shallow  depths, 
surface  effects  will  complicate  tha  aforementioned  sequence 
of  events.  When  the  shock  wave  hits  the  surface,the  atmos¬ 
phere  cannot  supply  appreciable  resistance  by  compression. 
As  a  result,  a  reflected  wave  with  a  negative  pressure 
satisfying  the  zero-pressure  condition  at  the  surface  is 
formed  (figure  4.1).  Thus,  the  resultant  pressure  observed 
is  the  sum  of  the  direct  and  reflected  pressures.  Therefore, 
the  reflected  wave  arriving  at  point  A  will  create  a  sudden 
drop  of  the  pressure  to  a  smaller  value.  This  is  known  as 
tha  "cut-off"  phenomenon,  typical  of  free  surface  effects 
(figure  4.2). 


I - - - i 


Figure  4.1  Surface  Effect  on  a  Shock  Wave. 


D.  PRESSURE  DETERMINATION 

As  detailed  in  a  previous  section,  the  pressure  decay  at 
any  point  is  roughly  exponential  so  that  it  can  be  written: 

P*  <A)  *  Pn(A)  qW  (4.1) 


It  has  been  found  that  the  fundamental  descriptors  cf  an 
underwater  explosion  attack  are  the  charge  size  (equivalent 
weight  of  TNT)  and  the  charge  standoff  (shortest  distance 
between  charge  and  target).  Theoretical  developments  about 
the  spherical  wave  detailed  in  [Sef.  10]  have  shown  that  the 
peak  pressure  Pm  at  any  point  can  be  reasonably  approximated 
by: 


Pm  =  K. 


1  (W^/B?3 


(4.2) 


where  W  is  the  charge  size  in  pounds  of  TNT  and  R  is  the 
standoff  distance  in  feet. 

It  has  been  shown  as  well: 


V3  V3  a2 
9  =  k2W  (w  /R)  2 


(4.3) 


Ki,K2rAi  'A2  arG  eaiPiricaHy  determined  factors  that  depend 
on  the  type  of  explosives  used.  Their  values  for  several 
types  of  explosives  are  given  on  table  V. 


Figure  4.2  Cat-Off  Phenomenon. 


TABLE  7 

Explosion  Parameters 


IIBX-l 

Tf.'T 

PENT 

NUKE 

Pnax  Kj 

22347.6 

22305 

24539 

4.32X1C6 

A1 

1.144 

1.18 

1.194 

M3 

Dacay  Constant  !U 

.OS  6 

.050 

.052 

2.274 

Ao 

.247 

.185 

.257 

_  ■>'» 

E.  THE  EXPLOSION  IH  EPSA 

EPSA  features  two  different  way  of  describing  underwater 
explosions: 

-A  discrete  fora  in  which  the  user  inputs  the  pressure 
history  at  a  finite  number  of  times. 

-A  functional  form  that  uses  equations  4.2  and  4.3  .The 
program  requests  then  the  various  coefficients  and  parame¬ 
ters  describing  the  explosion. 


TIME  (MSEC) 


Figure  4.3  Incident  Pressure  Decay. 
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The  presence  of  free-surface  effects  can  also  be 
accounted  for  with  the  input  of  a  cut-off  tiae  after  which 
the  incident  pressure  is  set  to  zero  (figure  4.3>. 
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T.  AHALYS IS  AND  RB£0L1S 


A.  MODELS  STODIED 

In  order  to  compare  the  results  of  the  numerical  anal¬ 
ysis  with  the  experimental  data,  all  attention  was  focused 
on  the  Explosive  Power  Meter  (EPtl>  model  for  which  field 
tests  had  been  conducted.  The  EPM  model  is  a  stiffened 
cylinder  with  end  ca  ps  whose  dimensions  are  given  in  figure 
5.1  . 


Figure  5.1  Explosive  Power  Meter. 
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Rather  than  Modelling  the  end  caps  with  additional 
sheets  ,  it  was  found  to  be  more  efficient  to  take  into 
account  the  behavior  of  the  end  caps  using  two  rigid  end 
blocks  (figure  5.1).  The  effects  of  the  explosion  on  the 
deformation  of  the  end  blocks  is  negligible  compared  to  the 
deformation  of  points  located  outside  of  the  and  blocks. 
Therefore,  it  was  assumed  that  the  lotions  of  the  end  blocks 
were  pure  rigid  body  displacements. 

In  order  to  gain  some  insight  into  the  influence  exerted 
by  the  stiffeners  and  the  end  blocks,  a  preliminary  analysis 
was  conducted  on  a  ring  stiffened  cylinder  without  end 
blocks  as  well  as  on  an  unstiffenad  cylinder  . 

In  addition  to  the  study  of  the  EPM  model,  the  influence 
than  the  location  of  the  stiffeners  had  on  the  deformations 
throughout  the  shell  was  evaluated.  By  performing  a  compara¬ 
tive  analysis  of  the  deformations,  it  was  intended  to  mini¬ 
mize  and  control  the  damage  caused  to  the  structure. 

The  cylinders  tested  were  subjected  to  an  explosion 
occurring  at  the  distance  R=  200"  from  the  cylinder.  The 
location  of  the  explosion  was  symmetrical  with  respect  to 
the  longitudinal  and  transverse  axis  of  the  cylinder  (figure 
5.2)  .  A  spherical  type  TNT  explosion  of  intensity  H=50  lb 
was  selected  for  the  study.  It  was  therefore  determined  by 
the  following  parameters  (chapter  17)  : 

Ax  =  1.  18  A2  =  -.185 
Kx  *  22505  K2  *  . 058 

The  symmetry  with  respect  to  the  xy  and  yz  plane  has 
been  taken  advantage  of  by  modelling  one  quarter  of  the 
model.  After  a  certain  amount  of  sensitivity  analysis  was 
performed  on  simple  grids,  a  finite  element  grid  consisting 
of  1517  nodes  and  1440  elements  was  selected  (figure  5.3)  . 
For  each  of  the  cylinders  studied,  the  time  step  chosen  was 
equal  to  ^t  s  3.1  O’*  s  .  The  explosion  process  was  studied 
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Figure  5.2  Explosion  Location 


FEB  Discretization  1517  nodes,  1440 


over  a  time  interval  of  800  time  steps,  that  allowed  for  a 
shock  and  after-shock  analysis. 

The  displacements  computed  by  EPSA  consist  of  a  combina¬ 
tion  of  rigid  body  displacements  and  pure  deformations.  The 
deformation  modes  are  of  significant  interest  since  they  may 
induce  buckling  and  even  lead  to  the  destruction  of  the 
structure.  As  mentioned  earlier,  the  very  stiff  end  blocks 
have  pure  rigid  body  displacements.  The  displacement  of 
each  node  of  the  end  block  was  subtracted  from  the  displace¬ 
ment  of  each  node  cf  the  corresponding  row,  giving  the 
component  of  the  displacement  due  to  pure  deformation. 

For  each  of  the  aforementioned  cylinders,  the  deformed 
geometry  and  the  color-filled  contour  plots  of  von  Mises 
stresses  as  well  as  normal  displacement  were  displayed. 
Using  identical  color  assignments,  the  progressive  gross 
evolution  of  the  previous  quantities  were  evaluated  so  as  to 
allow  for  a  comparative  analysis  of  the  evolution  of  phys¬ 
ical  parameters  throughout  the  shell. 

Color-filled  contour  plots  allow  for  a  global  representation 
of  a  quantity  distribution  and  have  been  found  extremely 
valuable  in  the  interpretation  of  the  results. 

For  printing  and  processing  reasons,  it  was  not  possible 
to  include  color  pictures  in  this  document.  Instead,  the 
contour  plots  of  the  physical  quantities  under  study  have 
been  included. 

B.  AHALISIS  OP  HIHG  STIFFS HBD  CYLIKDEBS  WITH  EHD  BLOCKS 
1.  S£H  Wodej, 

The  contour  plots  of  von  Sises  stresses  at  time 
steps  20,  60,  100,  140,  180  ,  200  are  provided  on  figure  A.1 
to  A. 6  .  As  the  shock  wave  hits  the  structure,  it  appears 
as  if  the  stresses  propagate  through  the  shell  and  reach 
their  maximum  value  fairly  quickly  in  50  time  steps.  After 


100  time  steps,  when  the  structure  is  fully  enveloped  by  the 
shock  wave,  the  region  close  to  the  end  block  becomes 
heavily  stressed  (figure  A. 4).  In  addition,  some  concentra¬ 
tion  of  stresses  at  the  locations  of  the  stiffeners  can  be 
seen  (figure  A.  5).  At  later  time  steps,  the  pressure  becomes 
decreased  and  there  is  a  relaxation  of  stresses.  However, 
the  region  close  to  the  end  block  as  well  as  some  spots 
located  around  tLo  stiffeners  remain  heavily  stressed,  which 
may  indicate  local  buckling  (figure  A. 6) . 

2.  Controlled  Deformations 

In  order  to  obtain  a  more  uniform  distribution  of 
displacements,  the  stiffeners  have  been  shifted  towards  the 
end  blocks.  The  time  history  evolution  of  the  displacements 
was  studied  over  an  interval  of  833  time  steps  for  the  EPM 
model  as  well  as  for  the  model  with  shifted  stiffeners 
called  EPM2.  The  ccntour  plots  of  normal  displacements  at 
time  steps  200,  400,  600,  800  for  both  models  are  provided 
on  figure  A. 7  to  A.  14 

The  EPS  model  shows  a  significant  concentration  of 
deformations  occuring,  even  at  late  time  steps  (figure 
A. 10),  indicating  a  possibility  of  buckling.  Although  unex¬ 
pected,  the  fact  that  the  region  close  to  the  end  blocks 
undergoes  the  most  severe  deformations  has  been  confirmed  by 
experimental  data.  A  possible  explanation  to  this  phenomenon 
is  that  the  inertia  of  the  cross-section  of  the  cylinder  is 
relatively  uniform  along  the  cylinder,  except  at  the  end 
blocks  where  it  jumps  to  a  much  higher  value.  This  disconti¬ 
nuity  in  the  inertia  results  in  concentrations  of  stresses 
that  cause  the  aforementioned  phenomenon. 

The  cross-section  inertia  of  the  EP(I2  model 
increases  more  smoothly  because  of  the  distribution  of  stif¬ 
feners  along  its  axis.  Thus,  the  concentration  of  stresses 
has  a  lower  magnitude  and  the  regicn  close  to  the  end  block 
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suffers  less  damage  than  in  the  EPM  case.  It  can  also  be 
seen  that  the  deformations  are  aore  uniforaly  distributed 
alono  the  axis  of  the  cylinder.  Above  all,  the  deforaations 
at  the  late  tiae  steps  are  not  as  large,  indicating  that  the 
chances  of  buckling  are  significantly  lower  for  the  EP32 
aodel  (figure  A.  14). 

Therefore,  by  perforaing  an  optimization  of  the 
location  of  the  stiffeners,  the  designer  can  counteract  the 
concentrations  of  stresses  and  the  buckling  phenomena  that 
occur  in  the  region  close  to  the  end  blocks.  It  is  believed 
that  controlled  defcraations  can  have  a  very  significant 
influence  on  the  survivability  of  a  structure  submitted  to  a 
shock  wave. 

C.  ANALYSIS  OF  ONSTIFFENED  AND  HINS-STIPFENED  CYLINDER 

It  was  decided  to  study  the  progressive  gross  responses 
of  an  unstiffened  cylinder  as  well  as  that  of  a  ring- 
stiffened  cylinder  without  end-blocks.  Both  cylinders  have 
the  same  external  dimensions  as  the  EPS  model.  The  ring- 
stiffened  cylinder  is  similar  to  the  EPfl  model  except  for 
the  fact  that  the  end-block  has  bean  replaced  by  a  standard 
stiffener.  The  evolution  of  von  flises  stresses  at  tiae 
steps  40,  80,  100  is  provided  in  figures  A. 15  through  A. 17 
for  the  unstiffened  cylinder.  The  evolution  of  von  Mises 
stresses  at  tiae  steps  40,  80,  100,  150  is  provided  in 
figures  A. 18  through  A.  21  for  the  ring-stiffened  cylinder 
without  end-blocks. 

For  the  unstiffened  cylinder  it  is  observed  that  the 
stresses  propagate  quickly  throughout  the  shell  and  that 
within  a  hundred  tiae  steps  an  instability  phenomenon  occurs 
showing  the  existence  of  local  buckling  (figure  A. 17).  At 
later  tiae  steps  the  buckling  spreads  over  the  entire  stiff¬ 
ened  shell. 
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The  evolution  of  von  fiises  stresses  in  figures  A.  13 
through  A. 21  shows  that  the  ring-stiffened  cylinder  can 
withstand  much  higher  stress  levels  than  the  unstiffened 
shall.  This  result  was  expected  since  the  stiffeners  provide 
the  stiffness  required  to  withstand  higher  loads.  At  time 
step  100  the  unstiffened  shell  is  already  subjected  tc  local 
instability  characterized  by  a  "hard  spot"  in  the  Biddle  of 
the  model  (figure  A.  17).  On  the  other  hand,  the  stresses  in 
the  stiffened  cylinder  are  much  more  evenly  distributed 
throughout  the  aodel,  with  high  amounts  of  stresses  concen¬ 
trated  around  the  locations  of  the  stiffeners  (figure  A. 21). 

It  can  also  be  observed  that  a  significant  concentration 
of  stresses  occurs  at  the  extremities  of  the  stiffened  shell 
(figure  •&.21).  Recalling  that  the  end-block  has  been 
replaced  by  a  standard  stiffener,  the  cross-section  of  the 
shell  has  a  greater  inertia  at  the  extremities  due  tc  the 
fact  that  the  two  stiffeners  located  at  the  extremities  are 
close  to  each  other.  Therefore  this  phenomenon  is  similar  to 
the  one  encountered  when  studying  the  EPM  model.  However 
the  concentration  of  stresses  for  the  ring-stiffened  cylin¬ 
drical  shell  is  less  significant  than  for  the  EPM  model,  due 
to  a  smaller  discontinuity  in  cross- section  inertia. 
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VI.  CONCLUSION 


— ' - — 


A  FORTRAN  module  that  merges  the  finite  element  code 
EPSA  with  the  color  graphics  system  PATRAN-G  has  beer, 
designed  and  succesfully  completed.  The  non-linear  elasto- 
plastic  responses  of  various  types  of  submerged  cylindrical 
shells  have  been  evaluated  using  the  EPSA/PATRAN-G  system. 

The  analysis  of  the  progressive  gross  responses  of  a 
ring- stiff ened  cylindrical  shell  with  end-blocJcs  is  believed 
to  provide  useful  information  regarding  the  behavior  of  a 
submerged  structure  subjected  to  an  underwater  explosion. 
The  influence  of  the  location  of  the  stiffeners  on  the 
deformations  has  been  studied  and  is  also  believed  to  be  of 
significant  help  in  the  determination  of  an  optimal  design 
that  will  minimize  the  damage  due  to  an  underwater 
explosion. 

The  utilization  of  the  color  graphics  system  in  the 
interpretation  of  the  results  of  analysis  has  been  found  to 
be  an  extremely  valuable  tool.  It  is  the  author's  belief 
that  the  use  of  color  graphics  systems  will  become  increas¬ 
ingly  important  in  the  analysis  of  complex  phenomena  such  as 
underwater  explosions  on  submerged  structures. 
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CPU  208  STEPS 
6  0000013E-06 


EPH  200  STEPS 
t  7'»‘>99$>2E-04 


200  STCP8 


PI1  200  STEPS 
6.3999969E-04 


26608 


200  STCPS 
000001 4E  04 


008  STEPS 
00600HC-04 


800  STEPS 
2000047E-03 


800  STEPS 
8000093E-03 


EPH  860  STEPS 
2 . 40081 39E-03 


800  STEPS,  STIFFENERS  SHIFTED  TO  END  BLOCK 
00000ME-04 


Sea  STEPS.  STIFFENERS  SHIFTED  TO  END  BLOCK 
200004FE-03 


8eeee93c-03 


880  STEPS,  STIFFENERS  SHIFTED  TO  ENO  BLOCK 
40081 39E-03 


Figaro  A. 17  Ton  Hises  stresses,  step  100,  unstif.  shell. 
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STIFFENED  CYLINDER,  NO  END  BLOCKS 
2eeeeB3E-B4 


>  STIFFENED  CYLINOER,  NO  END  BLOCKS 
3959979  C-0-4 


STIFFENED  CYLINDER,  NO  END  BLOCKS 
S>9$9S>6?E-04 


APPENDIX  B 

RE  VI  EH  OP  NONLI HEAR  FINITE  ELEMENTS 

A.  IITHODOCTION 

This  appendix  is  intended  no  give  the  reader  some 
insights  into  nonlinear  finite  elements.  The  reader  is 
assumed  to  have  some  previous  knowledge  of  finite  element 
theory.  The  basic  principles  of  the  theory  will  be  quickly 
reviewed,  but  the  study  will  focus  on  the  problems  that 
occur  when  dealing  with  nonlinear  theory.  Most  of  the 
information  has  been  takan  from  ^Ref.  6]  as  well  as  from 
the  course  the  author  had  at  M.I.T.  with  K.J.  Bathe  in  1982. 

B.  THE  NEED  FOB  A  N2H  THEORY 

Considering  a  coordinate  frame  defined  by  (i,  j,k),  a 
body  of  volume  V  in  which  point  A(xlrx2,x3)  is  subjected  to 
the  displacements  (ulru2  ,<13),  corresponding  to  a  strain 
vector  (e}  (figure  B.1). 

In  the  following  sections,  the  superscripts  0  and  t  will 
refer  to  the  body  at  time  0  and  t  respectively,  the 
subscripts  0  and  t  will  refer  to  the  configuration  at  time  0 
and  t  respectively.  This  chapter,  for  the  purpose  of 
simplicity,  will  first  deal  with  the  static  nonlinear  anal¬ 
ysis  of  the  material. 

In  the  linear  theory  of  finite  elements,  one  uses  the 
well  known  Cauchy  stress  tensor  associated  with  the 

engineering  strain  tensor  e^  ♦  3u^j  . 

Then,  the  principle  of  virtual  work  is  written: 

where  fcR  represents  the  virtual  work  of  externally  applied 

forces  and  <$e  is  a  virtual  strain. 


Figure  B.1  Geometric  Conventions 


Equation  (B.1)  is  then  discretized  over  the  body  and  becomes 
a  set  of  integrations  over  each  of  the  finite  elements.  In 
the  case  of  a  large  displacement,  the  volume  of  the  body 
over  which  the  integration  is  performed  might  have  signifi¬ 
cantly  changed.  Also  notice  that  equation  (B.1)  is  written 
in  rhe  original  coordinate  frame  (defined  at  t=3)  and  that 
and  ^  refer  to  the  current  configuration  of  the  body. 
The  Cauchy  stresses  at  time  t+ a  t  cannot  be  obtained  by 
adding  an  increment  due  to  the  straining  of  the  material  to 
the  stresses  at  time  t  .  The  rigid  body  rotation  of  the 
material  has  to  be  taken  into  account  since  the  Cauchy 
stresses  vary  under  rigid  body  motions.  Therefore,  we  must 
perform  the  integration  of  equation  (B.  1)  over  the  unknown 
current  volume  with  respect  to  the  original  geometry  that 
could  be  significantly  different  from  the  current  one. 

The  above  discussion  emphasizes  the  need  for  a  new  set 
of  stress  and  strain  tensors  that  would  alleviate  the  afore¬ 
mentioned  problems  and  enable  the  integration  of  the  prin¬ 
ciple  of  virtual  work  to  be  performed. 


C.  DEFINING  NEW  STBESS  AND  STBAIN  TEMSOBS 


1.  Green-Laqrange  Strain  Tensor 

The  structure  introduced  earlier  is  considered  and 

two  points  A  and  B  at  t=0,  of  coordinates°A  (°x)  ,  °3(0x) 

are  defined  .  At  time  t,  the  body  has  been  deformed  and  A 

and  B  have  moved  to  tA(tr)  and  tB{tr)  (figure  B.2). 

(  ( 


Figure  B.2  Displacements  Conventions. 


A  Taylor  expansion  is  used  to  express  the  coordinates  of  B 
as  a  function  of  the  coordinates  of  A. 

txi  =  ^  ♦  9  txi  (°xj-  °tj  (B.2) 

3  *j 

or  with  d  txi  -  txi  and  d°xi*  °xi  -  °*i  ,  gives  : 

(d  ^  =  <3^  )  (d  Sci)  (B.3) 

3  Xj 

{ d fcx}  =  [  qX  ]  (d°x}  (B.«) 


where  [  o*  ]  3  ( 

dX) 

cd^j 

3  (3  » 

(d°x|  =  (d°Xi) 

In  other  words. 

aquation 

(B.4) 

expresses 

how  a  small  fiber 

defined  by  the  vector 

{d  °X) 

at  t=0 

has  rotated  and 
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extended  between  time  0  and  time  t  when  it  becomes  {dtx} 
The  matrix  [qX]  is  called  the  "deformation  gradient"  . 

The  new  length  fcdS  of  the  fiber  will  be  : 


, ,  t  ,  T  , 

(d  x}  {d  x} 


<  dS)*  = 
matrix) 

and  therefore 

( ^ds)  2  =  (C  ox  1  C^°x} )T  (Co  x  ]  Cd°x} ) 


(  T  refers  to  the  transposed 


(fcdS)2  =  (d°xj  CqC]  {d°x) 


(B.5) 


t  t  T  t 

Cocl  =  Coxl  Cox3  is  called  the  "Cauchy-Green  deformation 
tensor" . 

Notice  that  if  the  Cauchy-Green  deformation  tensor  is  iden¬ 
tity,  equation  (B.5)  indicates  that  the  length  of  any  fiber 
will  not  vary.  In  other  words,  whenever  rigid  body  motions 
are  considered,  the  Cauchy  Green  deformation  tensor  is  iden¬ 
tity  since  the  fibers  do  not  stretch. 

The  principles  of  the  finite  element  method  will  now 
be  quickly  recalled.  Assume  that  the  displacement  (i^  )  of 
any  point  of  the  body  can  be  written  as  a  combination  of  the 
displacements  of  N  selected  points  called  "nodes"  : 


k=l 


(B.6) 


Where  the  h  are  interpolation  functions  that  depend  only  on 
the  geometry  of  the  body.  In  addition,  the  nodes  are  chosen 
so  as  to  get  a  division  into  quadrilateral  elements  and  it 
is  assumed  that  the  displacement  of  any  point  is  only  a 
function  of  the  displacements  of  the  corner  nodes  of  the 
element  it  beLongs  to.  Then  equation  (B.6)  becomes 
4 

(Ui)  =  £  ^  (Ui)k  (B  .7) 

k=*l 

Becalling  that  (u^)  is  simply  (^i)  -  (°Xi)  gives: 
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(B.  8) 


( 


\i  -  <°\>  *1  V“i>* 


k=l 


The  deformation  matrix  [ Q  X]  can  be  expressed  very  simply 
using  the  previous  equation  : 


(B.9) 


Define  row  the  "Green  Lagrange"  strain  tensor  by  : 


Coe]  =  1/2  (Co  C  3  -  [i]) 


(B.  10) 


Sfhere  [I]  is  the  identity  matrix. 

Prom  the  previous  derivations,  it  can  be  observed  that  the 
Green  Lagrange  strain  tensor  is  0  for  rigid  body  motions. 
The  Green-Lagrange  strain  tensor  refers  to  the  body  at  time 
t  with  respect  to  the  initial  configuration.  This  is  why  it 
will  be  so  useful  in  dealing  with  large  deformations. 

Recall  that  the  ultimate  goal  is  to  apply  the  prin¬ 
ciple  of  virtual  work  to  the  structure  under  study.  In 
particular,  having  defined  a  new  strain  tensor,  the  relation 
giving  the  virtual  Green-Lagrange  strain  tensor  corre¬ 
sponding  to  a  virtual  displacement  (6  u|  must  be  known.  In 
the  case  of  a  linear  problem,  the  virtual  engineering  strain 
tensor  would  be: 

*  1/2  (  3uj +  3uj ) 

J  3xj  5“sq. 

In  the  case  of  large  displacements,  the  Green-Lagrange 
strain  tensor  should  be  used.  It  is  shown  in  [  Bef .  6]  tha+ 
the  virtual  Green  Lagrange  strain  tensor  corresponding  to 
virtual  engineering  strains  is: 


5  o*ij 

or: 


Sx_  Stem 


(B.  11) 
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<5  t®ij  *  3*i  3rj  $0£ij  (3.12) 

3xm  3xn 

Having  defined  a  strain  tensor  which  is  invariant  under 
rigid  body  notions,  a  stress  tensor  corresponding  to  the 
Green-Lagrange  strain  tensor  needs  to  be  defined. 

2.  Stress  Measures 

Starting  with  the  Cauchy  stress  tensor  ,  the 

Piola-Kirchof f  stress  tensor  is  defined  : 

O^ij  S^iui  tTmn?.xjn  (B.13) 

Where  %  ,  fcp  are  the  densities  .of  the  material  at  time  0  and 
t  respectively  and  (  °tXij  =  [  °X  ]l  is  the  inverse  of  the 
deformation  tensor  defined  previously. 

Equation  B. 13  can  be  easily  rewritten  in  equivalent  form: 

"'mn  Oxim  BSij^jn  (B.1U) 

It  can  also  be  shown  by  using  the  principle  of  mass  conser¬ 
vation  that  : 

Op  =  t0  det[  §X  ]  (B.  15) 


D.  PRINCIPLE  OP  VIRTUAL  WORK 


The  principle  of  virtual  work  of  equation  (B.1)  is 
rewritten  using  the  strain  and  stress  tensors  defined  previ¬ 
ously  in  equations  (B.12)  and  (B.14|  : 


t  R 


or: 


0  0*nk  £  *iin  5  eij  ^ 


(B.  16) 
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ur  =  nz 


I  Op  oS^i^ 


and  using  equation  (B.15) 


^  =  /  0 Si jd osi jd V 


(B.  17) 


(B.  18) 


Thus,  the  principle  of  virtual 
expressed  in  terms  of  a  new  set 
tensors,  integrated  over  the  original 


work  has 

of  stress 

o„ 

volume  V 


been 

and 


simply 

strain 


E.  THE  IHCRESE3TAL  CONTINUUM  MECHANICS  EQUATIONS 


In  this  section,  the  principle  of  virtual  work  will  be 
applied  to  the  structure  and  the  incremental  formulation 
using  the  Piola- Kirchcf f  and  Green- Lagrange  tensors  will  be 
developed  .  Non  linear  terms  will  arise  from  the  rather 
complicated  definition  of  the  strain  tensor,  but  if  must  be 
pointed  out  that  the  new  formulation  provide  the  means  for 
the  modelling  of  large  deformations. 

Assume  that  the  configuration  of  the  body  at  time  t  is 
known,  the  configuration  at  time  t  +  At  must  be  determined. 
Writing  the  principle  of  virtual  work  at  time  t+it  gives  : 


t+A{?  :  external  work  at  time  t  +  At 

e  Tat 

5  oGij  virtual  increment  in  G.L.  strain 
o^ij  :  stress  state  at  time  t+At 


(B.  19) 


Separating  between  the  known  terms  that  refer  to  the  config¬ 
uration  at  time  t  and  the  unknown  terms  that  are  the  incre¬ 
ments  of  stress  and  strain  between  time  t  and  time  t+&t 
gives  : 
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I 


t+  At  t 

0S  =  o  ^ij  +  o  Sj.j  ar.d  z 


0  eij  +  0  G  ii  G 


(B.  20) 


o Sij  andoGij  are  simply  the  increments  in  stress  and  strain 
respectively. 

It  is  derived  in  [Ref.  6]  that  the  increment  in 
Green-Lagrange  strain  o^ij  is  made  of  a  linear  part  oeijand  a 
nonlinear  one(fiij  .  The  term  linear  refers  to  the  increment 
in  displacement  u^  . 


o  Gij  =  o^j  ♦oOij  cr  5  oPij  =  5  (  oeg  +0TUj  )  (B.21) 


Osing  equation  (B.21)  in  equation  (B.  19)  gives  : 
f  (  $5%  *  (  5  (0eij  ♦o'lij  ))  dV  = 

Jn  _ 


.  (B.  22) 


Again,  separating  between  the  Known  and  unknown  terms  gives 


[  QSi^ijdV  QS  5QHij  -  ~y^0^ij^Oeij 


(B. 23)  . 


'i  i 

L 


For  small  increments  in  displacements,  equation  (B.23) 
written  at  time  t  indicates  that  5  0eij  =  Soey  .  This  signi¬ 
fies  that  in  equation  (B.21)  the  non  linear  term  is  negli¬ 
gible 

Then  the  constitutive  law  of  the  material  detailed  in 
chapter  II  allows  to  relate  stresses  and  strains: 


O^ij  ~  0*i.jnsC£ij  Cpijcs  0®ij 


(B.  24) 


and  B.23  becomes: 


jWj  cftjscasdV  ♦j'SVo^dV-^R  cPijdv 


(B.  25) 


» 


The  right  hand  side  terms  of  the  previous  equation  are 
known : 

is  the  linear  increment  in  virtual  strain  involving  only 
known  terms 

Bsij  is  known  from  the  previous  time  step 

tS-At 

R  is  the  work  of  external  prescribed  virtual  forces. 

The  left  hand  side  of  B.  25  is  unknown  since  0®i.jand  ^o'lj 
involve  the  displacement  from  time  t  to  time  t+At  . 

P.  PIUTE  ELEMENT  DISCRETIZATION 

Equation  (B. 25)  will  be  discretized  over  the  structure, 
using  the  finite  element  approximation  defined  in  section 
B.9  .  Let  N  be  the  number  of  nodes,  the  principle  of 

virtual  work  will  be  invoked  N  times,  setting  a  unit 
displacement  at  each  node  in  turn: 

5  u  k  =  1  ,  5Uj  *  0  ,  k#j 

A  system  of  equations -whose  unknowns  are  the  nodal  displace¬ 
ments  is  obtained  Let  "  (  A  u}  be  the  vector  of  unknown 
displacements,  {PJ  be  the  vector  of  nodal  point  forces 
equivalent  to  the  internal  stresses. 

Then  equation  (B.25)  can  be  rewritten  in  matrix  form  : 

CXHAu}  *  (B.  26) 

[  qK  L]  and  C  qkm,]  are  knovn  from  the  material  characteristics 
at  time  t  and  correspond  respectively  to  a  linear  and  nonli¬ 
near  contribution.  It  is  therefore  possible  to  solve  equa¬ 
tion  (B.25)  for  {Au}  .  However,  because  of  the  assumption  in 
equation  (B.24),  the  exact  solution  might  not  be  reached 
immediately.  Purthermore,  depending  on  the  time  step  size, 
the  solution  process  might  even  be  unstable!  In  any  case. 
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.an  iteration  for  solving  equation  (B.26)  must  be  performed 
until  the  exact  solution  is  reached  . 

A  widely  used  scheme  is  the  "modified  Newton  iteration" 
defined  by  the  following  aquation  and  boundary  conditions  : 


(  CoKl]  ♦  CoKniJ  )  (Au)i 


=  r 1 


(B.  27) 


and 


=  {* 


(Au£ 

=  ^u)  and  {  ] 


with  the  initial  conditions  : 


and  {  F}°  =  £?} 


iAU}1  is  the  vector  of  incremental  nodal  point  displacements 
at  iteration  i. 

is  the  vector  of  applied  loads  (constant  in  the  itera¬ 


tion) 
tt-A:  i-l 


{  F}  is  the  vector  of  nodal  point  forces  equivalent  to  the 
stresses  at  time  t>it,  iteration  i-1  . 

At  the  first  iteration,  equation  (B.27)  reduces  to  equa¬ 
tion  (B.26)  giving  an  increment  of  displacement  (Au)1.  Then 
a  better  approximation  of  oers  *  5 cr±i  is  obtained  .  gSri  is 

J  fclAt-  J 

updated  to  the  new  state  of  stresses  and  becomes  o6ij 
Equation  (B.27)  is  then  used  to  determine  the  new  increment 

in  displacement  (Au)2  ,  and  so  on  unt.ll  the  increment  in 

t a-  &  t-fct  i 

{  8}  *  C  F  }  in 


displacement  is  small  enough,  so  that 
equation  (B.27) 


G.  INCLUSIOB  OP  DIN  ABIC  PORCES 

If  the  loads  are  applied  rapidly,  inertia  forces  need  to 
be  considered  and  a  truly  dynamic  problem  has  to  be  solved. 
Using  d'Alembert's  principle,  the  eLement  inertia  forces  are 
simply  included  as  part  of  the  body  forces.  Let  (u)  be  the 
vector  of  nodal  accelerations  and  [  S  ]  be  the  mass  matrix  of 
the  system.  Then  the  principle  of  virtual  work  is  written  in 
the  following  way  : 
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t  t 

(Co^Cau}  ♦  CqK^Cau  }  *  B  -CMJCau}  -f  (B.28) 

Equation  (B.28)  represents  a  system  of  differential  equa¬ 
tions  of  second  order.  If  the  non  linear  term  (J^kKL  ]  were 
negligible,  the  solution  could  be  obtained  by  standard 
procedures  for  solving  differential  aquations  with  constant 
coefficients.  However,  the  procedures  proposed  for  the  solu¬ 
tion  of  general  systems  of  differential  equations  can  become 
very  expensive  if  the  order  of  the  matrices  is  large. 
Therefore,  whenever  the  system  is  linear  or  nonlinear,  seme 
effective  methods  for  solving  the  equations  governing  the 
equilibrium  are  required. 

1 .  Direct  Integration  Methods 

The  essence  of  direct  integration  methods  is  based 
on  two  ideas.  First,  it  is  aimed  to  satisfy  B.28  only  at 
certain  time  intervals  apart.  Second,  a  variation  of  accel¬ 
eration  velocities  and  displacements  is  assumed  within  each 
time  step.  The  form  of  the  assumption  determines  the  accu¬ 
racy,  stability  and  cost  of  each  method. 

In  the  following,  assume  that  the  initial  conditions 
(displacements,  accelerations,  velocities)  at  time  0, 
denoted  (  °u  ,°u  ,°  ii  )  are  known.  In  the  solution,  the  time 
span  under  consideration,  T  ,  is  subdivided  into  n  equal 
time  intervals  At.  Assuming  that  the  solution  is  known  at 
time  t  ,  the  methods  of  getting  the  solution  at  time  t+At 
will  be  investigated. 

2.  Central  Difference  Method 

In  the  Central  Difference  method,  a  finite  differ¬ 
ence  approximation  will  give  the  acceleration  at  time  t  : 


i  Js-At  _t  tsAt  . 

-  2  u  ♦  u  ) 


(B.  29) 
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iriting  the  principle  of  virtual  work  at  time  t  and 
substituting  into  equation  (B.30)  gives  then: 


is -tk. 

i^]{u} 


♦  CoK  ]Cfcu}  =  {tB} 

=  {tH}  -  (C<? K  ]  -  2[  S  ]/A  t2  )  {%} 


tJt 

a]{  u} 


(B.  30) 


(B.31) 


The  previous  equation  gives  the  da  formation  at  time  t+At 
from  the  characteristics  of  the  system  at  tiae  t  . 

When  [M]  is  diagonal,  which  is  frequently  the  case 
for  mass  matrices,  the  solution  at  time  t+  4 1  does  not 
involve  any  triangular  factorization  of  the  matrix  [M]  , 
thus  leading  to  mor9  efficient  computations. 

The  shortcoming  in  the  use  of  the  central  difference 
method  lies  in  the  time  step  restriction:  for  stability,  the 
tiae  step  size  t  must  be  smaller  than  a  critical  time  step 
#crwhicL  is  equal  to  I  4  ,  where  T  is  the  smallest  period 
of  the  finite  element  system. 

The  central  difference  scheme  is  fairly  easy  to 
implement  for  the  integration  of  a  system  of  non  linear 
differential  equat ions.  However,  because  of  the  limitations 
of  the  time  step  ,  it  might  not  be  suitable  for  cases  when 
loads  are  varying  at  a  slow  pace. 

3.  Implicit  Integration  Schemes 

Since  the  interest  of  this  study  lies  in  central 
difference  schemes,  this  section  will  be  limited  to  a  short 
description  of  the  fundamentals  of  implicit  time  integration 
schemes. 
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■V 


Implicit  time  integration  schemes  use  the  principle 
of  virtual  work  written  at  time  t+  4t  and  not  at  time  t  as 
for  the  central  difference  method  : 

Hit  t-ist  t^t  HAfc 

[  M]  (ii)  ♦  CoK]  {u}  =  (  H}  (3.32) 

Hit 

Again,  using  a  finite  difference  approximation  of  (u)  and 

tH^t 

replacing  into  equation  (B.32)  enables  to  solve  for  {  u)  . 

bhit 

Since  the  formulation  involves  the  rigidity  matrix  fo  X  ]  and 
tha  erternal  work  {  R)  which  are  both  unknown,  the  system 
has  to  be  solved  in  a  similar  way  to  the  Modified  Newton 
iteration  that  was  detailed  previously. 

Implicit  time  integration  schemes  are  stable  regard¬ 
less  of  the  size  of  the  time  step  used.  However,  if  the  time 
step  size  is  too  large,  significant  errors  can  be  accumu¬ 
lated  at  each  time  step,  leading  to  unrealistic  results. 

The  reader  will  find  more  details  on  the  various 
implicit  method  in  [  Hef .  6].  Yet,  it  can  pointed  out  that 
implicit  methods  are  more  tedious  to  implement,  on  the  other 
hand,  a  larger  time  step  can  be  used  in  the  solution  proce¬ 
dure,  which  can  be  of  extreme  importance  when  studying 
phenomena  over  a  significant  period  of  time. 
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IP  PEN PIT  C 

HOW  TO  (JSE  THE  TRANSLATOR  EPS  A-P  ATE  AN-G 


This  appendix  is  intended  to  explain  in  denail  the  use 
of  E PS A/P AT RAN -G  post-procsssing  facilities.  It  is  divided 
in  three  sections  that  will  deal  successively  with  1)  the 
displaying  of  the  original  model  ;  2)  the  deformed  geometry; 

3)  the  contour  plots  of  element  and  nodal  points  quantities. 

A.  DISPLAYING  THE  OBIGINAL  MODEL 

When  making  an  initial  EPSA  analysis  on  a  particular 
structure,  the  geometry  of  the  model  has  to  ba  input  into 
PATRAN-G.  As  explained  in  chapter  III,  all  the  geometrical 
information  is  contained  in  a  file  FOS019.DAT  that  is 
created  each  time  an  EPSA  run  is  aade.  The  input  of  the 
original  geometry  must  be  made  via  the  neutral  input  mode  of 
EPSA.  The  procedure,  starting  from  the  "logon"  to  PATRAN-G 
is  the  following  ; 

-  Select  the  GO  option 

-  Select  the  new  data  file  option  (option  1) 

-  Select  the  neutral  system  (option  4) 

-  Select  the  input  mode  (option  4) 

-  Input  the  neutal  file  name  ;  FOR019.DAT 

The  original  geometry  will  then  be  displayed  on  the  screen 
It  is  often  found  convenient  to  have  a  perpective  view 
of  the  model  under  study.  In  than  case,  the  user  should  ; 

-  Issue  the  VIEW  command 

-  Select  the  rotation  about  the  absolute  axes  (option  1) 

-  Input  an  angle  of  rotation 

(23,-34,0  will  give  a  very  nice  view  but  any  angle  can  be 
input) 
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-  Issue  a  PLOT  command  to  have  the  model  displayed  ir.  the 
new  a  xe  s . 

An  example  of  the  procedure  is  provided  on  table  VI. 


TABLE  71 

Pinite  Element  Model  Input  Procedure 


MODE'?  1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3. DISPLAY  4. NEUTRAL  SYS.  5. END 

>4 

NEUTRAL  FILE?  i. CREATE  OUTPUT  2. INPUT  MODEL  3. POST-PROCESSING  4. END 
>2 

INPUT  NEUTRAL  FILE  NAME 
/FOR019.DAT 

DO  YOU  WISH  TO  OFFSET  ANY  NEUTRAL  INPUT  IDS?  (Y/N) 

>N 

EPM  200  STEPS  .NO  STIFFENERS  U*30. 

SHALL  UE  PROCEED  UITH  THE  READING  OF  THIS  FILE?  (Y/N) 

> 


i _ i 


The  PLOT  command  can  be  issued  anytime  to  display  the  orig¬ 
inal  geometry  on  the  screen. 

When  studying  complex  models,  one  does  not  want  the 
element  and  node  numbers  to  be  printed  along  with  the  geom¬ 
etry  of  the  structure.  The  command  SET,  LABS3,  OFF  followed 
by  a  PLOT  will  display  the  original  geometry  without  any 
labels  printed. 

When  a  model  has  been  input  into  PATRAN-G,  a  data  file 
PATRAN.DAT  is  created  on  the  user's  directory.  When 
connecting  with  PATRAN-G  at  a  later  time,  the  user  can 


90 


select  the  option  "last  data  file  "  (option3)  to  have  the 
original  geometry  displayed  on  the  screen  without  having  to 
input  the  model  again. 

B.  DEFORHED  GEOHETRY 

On  the  second  card  of  the  PATRAN-G  input  deck,  the  user 
specifies  the  number  of  PATRAN-G  displacement  data  file  and 
the  number  of  element  results  data  file.  Two  non-zero  inte¬ 
gers  in  free  format  must  be  placed  at  the  end  of  the  second 
card  (section  II  in  the  user's  manual)  reguesting  the  number 
of  displacements  and  of  elements  files  respectively. 

Assuming  that  the  user  has  made  a  200  steps  run  with  10 
output  requests  for  PATRAN-G  displacement  files,  ten  (10) 
files  FOR0 1 8. DAT  will  then  be  created  at  equal  time  inter¬ 
vals.  The  deformed  geometry  corresponding  to  time  step  100 
will  therefore  be  contained  in  FOR01 8 . DAT; 5. 

To  display  the  deformed  geometry  corresponding  to  time 
step  100,  the  user  should  issue  the  following  comm?-ds: 

-  RUN  ,DEF  ;  requests  deformed  geometry  option 

-  Input  the  name  of  the  displacements  file  :FOR018.DAT  5 

-  Select  the  PLOT  option  (option  3) 

-  Select  the  undefcrmed  geometry  (2)  followed  by  the 

deformed  geometry  (3).  An  example  of  the  procedure  is 

provided  on  table  711.  The  undeformed  geometry  superposed 

with  the  deformed  geometry  will  taen  be  displayed  on  the 
screen. 

C.  POST-PROCESS I«G  OF  INAL YIS  RESULTS. 

Element-related  quantities  liks  von  Jlises  stresses  are 
contained  in  P0R016.DAT  files,  nodal  point  quantities  are 
stored  in  FOR018.DAT  files.  As  described  in  chapter  III, 
each  column  of  those  files  contains  a  specific  quantity 


TABLE  VII 


Deforced  Geometry  Procedure 

> 

MODE1?  1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3. DISPLAY  4. NEUTRAL  SYS.  5. END 
>RUN, CON,COL,4 

CURRENT  FILE  FOR  NODAL  RESULTS  IS  FOR018.DAT 
NEUTRAL  RESULTS  FILE?  l.NEU  FILE  2. CURRENT  FILE 
>1 

INPUT  THE  RESULTS  FILE  NAME: 

>FOR0 18 . DAT ; S 
DATA  UIDTH  *  S 
FILE  TITLE  *EPN  800  STEPS 

7.5000129E-04 


DATA  VALUES  RANGE  FROM  -0.S2SE+00  TO  0.134E+00 
ASSIGNMENT*?  l.AUTO  2. MANUAL  3. SEMI-AUTO  4. USE  CURRENT  LEVELS  5. END 

>4 

PREVIOUS  CONTOUR  LEVELS  USED. 

MODE?  1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3. DISPLAY  4. NEUTRAL  SYS.  5. END 
>RUN,HI,C 


(i.e.  column  31  -of  FOH016.DAT  contains  the  von  Mises 
stresses).  The  reader  -will  refer  to  chaprer  III  for  the 
detailed  organization  of  those  files. 

Again,  assume  a  200  time  steps  analysis,  with  10  output 
requests  for  element  results  files.  The  user  might  want  to 
display  the  contour  plots  of  von  Mises  stresses  at  time  step 
100  .  In  this  case  the  following  commands  should  be  input  : 

-  RON,  CON,  COL,  31:  tells  PATRAN-G  to  look  at  column  31 
that  contains  the  von  Mises  stresses. 

-  Input  the  file  name  FOR016.DAT  5 

-  PATRAN-G  will  then  ask  for  a  color  assignment  (automatic, 
manual,  semi-automatic,  current  levels  used)  that  the  user 
will  select  according  to  his  needs. 

The  contour  plots  are  then  ready  to  be  displayed  : 

-  RON,  HI,  FR  would  display  the  color-filled  contour  plots 

-  RON,  HI,  CON  would  display  the  contour  plots  (color  lines) 
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The  von  Mises  stresses  are  the  lost  useful  element  quan¬ 
tities  to  be  displayed,  but  other  element-related  quantities 
detailed  in  chapter  III  like  the  x  and  y  stresses  in  local 
or  global  coordinates  could  be  displayed  as  well  by  looking 
at  their  corresponding  column  in  the  element  results  data 
file. 

Dealing  with  nodal  point  quantities,  the  displacement 
normal  and  the  velocity  normal  to  the  shell  are  very  mean¬ 
ingful  quantities  in  an  analysis.  They  are  respectively 
stored  in  column  4  and  5  of  the  displacement  results  files 
FORO  1  8.  DAT  . 

A  contour  plot  of  the  normal  displacement  at  time  step  100 
would  then  be  obtained  via  the  following  commands: 

-  RON,  CON,  COL,  4  (look  at  column  4) 

-  FOR018.DAT  5  (name  of  the  file* 

-  Color  assignment  chosen 

-  RON,  HI,  C  or  RON,  HI ,  FR 

Notice  When  the  fluid-structure  interaction  is  ON,  the 
normal  displacement  contained  in  column  4  corresponds  to 
pure  deformations.  the  rigid  body  contribution  having  been 
taken  out. 

All  the  element  results  processing  is  implemented  in  the 
routine  NEOSTRE,  all  the  nodal  points  processing  is  imple¬ 
mented  in  routine  NEODISP.  It  shouli  be  pointed  out  that  any 
modification  to  the  capabilities  of  the  translator  (i.a. 
being  able  to  display  other  types  of  quantities)  can  be 
made  by  modifying  these  routines  only. 
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ip  pen pi i  g 

LISTINGS 


This  appendix  contains  the  various  files  that  constitute 
the  translator  module.  The  submodule  that  displays  the  orig¬ 
inal  geometry  is  listed  on  the  first  four  pages.  It  is 
imbedded  in  file  FRANK. FOR.  The  submodule  that  takes  care  of 
the  post-processing  facilities  is  imbedded  in  file  DISP.  FOR 
and  is  listed  in  the  remaining  pages.  It  has  been  mentioned 
previously  that  2PSA  had  been  slightly  modified  to  accomo¬ 
date  color  graphics  capabilities.  The  only  interaction  of 
EPSA  with  the  translator  occurs  via  subroutine  calls.  All 
the  "calls"  occur  in  COMPOTE  (for  post-processing)  and 
REPORT  (for  the  original  geometry)  .  A  labelled  COMMON 
called  FRANK  has  been  created  and  is  defined  in  the  routine 
AAA  as  requirad  by  EPSA.  The  requests  for  PATRAN-G  outputs 
are  echoed  in  the  EPSA  output  file,  all  the  modifications 
for  that  purpose  having  been  made  in  routine  BEADIN.  The 
user  must  be  awara  that  the  size  of  blank  COMMON  array  A  has 
been  increased  to  store  the  deflections  in  x,  y,  z  direc¬ 
tions  instead  of  only  the  z  direction  previously. 
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subroutine  Orel  in 
dimension  4(1) 
coi»o"  i  a  ( 1  ) 

eauivalence  (ta(l).a(l)) 

connon/ssite/  ibgCll/nqj.neltot.nlbd.nload.nbrect. 

noauadr  isheet.  norofs.  nsots.  nstiDts.  nvots.  nhots. 

1  nsjtvo<  nn j ,  nntot.  Iqdso.  liaud»  lbcatc 
cannon  /cotrj/  nsteos.  noeg.  nend>  nsheet 

cannon  /stab  /  i b t  (  1  )  . i s s i t e »  i  soa r  .  >  v e I o  <  i s t re » j * n a s » i i e 1 m , j bn  a t 
jllod»jloco.joret»jlhis»isrrn,iforc»j*loc»jnqi,jnni,jnaceg. 

1  j 1 s i de 

connon/frank/nfntot  .  1  1 u 
COMMON/ T ITRE/NT I TLE (80) 

CHARACTER*”  RUFF 
CHARACTER'S  TITLE 
CHARACTER'S  TIM 

character *r  ver 

this  routine  comoutes  new  nu«ber  of  nodal  ooints 
for  anv  given  sheet 

kJ=JNGBEG-JNNI 

nfntot=nntot-U( JNNI )-IA( JNQBEG-1 )-2«(KJ-2) 

I  1 t  voe  =  2S 
11 kc=l 
1  1  i v  =  0 
1  1  i  0  =  0 
I  1 nlsO 
I ln2  =  0 
I  In3  =  0 
1  lnuso 
1 InSsO 

write(llu.lO)  Ittvoe.l lid. lliv. like. llnl,l)n2,llnS.11nU,lln5 
fornat(i2«8i8) 

writed lu.ll)  (((TITLE  (I).  I  si.  80) 
f arna  t ( 80  Al) 

taking  care  of  second  oacket 

1 1 t  yoe=2b 
1 1 kc  =  l 
1 lnlsnfntot 
1 1 n2=nel tot 

-ritedlu.10)  lltyoe.llid.lliv.llkc»llnl,lln2,llnS.llnU,lln5 

call  aateCBUFF) 

wr i t  e ( 1 1 u» 1 2)  BUFF 

fornat (»,3*. ’ 17: 12:09',5X, • 1 .«• ) 

return 

end 


Subroutine  shconn 
di mens  ion  A ( 1 ) 
connon  i a ( 1 ) 

eauivalence  (ia(l).a(l)) 

conmon/ssite/  iba(l).naj»neltOt»nlbd.nload.norect. 

1  noauadr  isheetr  norotSi  nsotSi  njtrotSi  nvots.  nhots. 

2  nsstvoi  nn j .  nntot.  Igrtso.  lidud.  Ibcelc 


qmQP 


AD-A139  072  UNDERWATER  SHOCK- INDUCED  RESPONSES  OF  SUBMERGED 

CYLINDRICAL  STRUCTURES(U)  NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY  CA  F  G  DAUBE  DEC  83  ' 

UNCLASSIFIED  F/G  20/4 


o  n  o 


\ 


c  common  /coira/  nsteDii  nneii  nend#  nsh*et 

C3**on  /stab  /  ibt(l)»jssije,jsoar,jvelo»isrre.jimas(jielm,jbmat. 

1  j11cd»ilo3o»iopet»i1His,j3trn,iforCfi»loCfjnqi,jnni,jnooeg» 

2  j lside 
C 

cornmon/frank/ntntor  <  1  lu 

e 

c  this  routine  will  write  the  elegant  corner  noties  ot  each 

c  element  on  neutral  tile 

c 

1 1 1  ype=2 
1 1 Vc  =  2 

1 1 noa=a 
1 1 iv=« 

LLN1 sO 
LLN2=0 
LLN JsO 
LLNasO 
LLN5=0 
LLC0NF=0 
LLPI0=0 
LLCE I =0 
THET1=0. 

THET2=0. 

THETJsO. 


INITIALIZATION  OONE 

norev=0 
1  1  e 1  *  j  nni - j  noi 

do  200  k*l  ,  1  lei 
1 1 rowsl A ( j nai tk-1 ) 
c  number  o*  elts  in  each  row 
do  100  j  s 1 , 1 1  row 
let  ®IA( jnqbegth-1 ) ♦ j-l 
n 1  el  1  * i rnorev 
n 1 e 1 2s i ♦ 1 *nprev 
nlel3sj»l+norev»l I  row* I 
nlel4=i*norev*l I row*l 
c  ready  to  disolay  oacket 
e 

1 1 id=LEL 

write!! 1 u#  80 )  1 1 tyoe, 11  id. 1  I i v , 1  I ke .LLN1 ,LLN2, LLN3, LLN« ,LLN5 

80  f orma t ( i  2 »  8  f  8  ) 

write(llu.8l)  1 lnod,LLC0NFrLLBI0,LLCEI, THET1 ,THET2, THET3 

81  f orma t ( i 8 > 3 i 8 , 3e 1 6 . 9 ) 

writed 1 u»  82 )  nlel I»n|»12,nlel3fnle14 

82  format(10i8) 

1 00  cont i nue 

e 

norevsnorev*llrow*l 
200  continue 

return 

end 


subroutine  sheetf 
dimension  All) 
common  i a ( 1  ) 

equivalence  Ci»Cl)»aCl)) 
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n  n  n 


eo»*on/jjij*/  ibq(!)*naj#neltot.nlbi,nloaa,nbreet. 

1  noouad#  isheet.  noratSi  nsots.  nstrots.  nvots.  nhpts. 

2  nsstva#  nnj,  nn tot.  Ijdso.  lia ud.  lbealc 

c  commo n  /eoara/  nsteos.  nbea.  nend,  nsneet 

co»«on  /srab  /  i b t ( 1 ) » i ss i te, j soar , j we  1 o< i s t re .  j <*as , j i e 1 m» j bma t , 

1  i I lori. i I odo . joret . j 1 h i s . j  st  rn. j  fore. j  x 1 oe . jnai . j  nnt . j  nooeg. 

2  i 1  side 

e 

common/frank/nfntot.  1  lu 
eharaeter*l  qtyoe 
e 

c  renumbers  the  nodes.  chanae  coordinates*  display  oacket  1 

1 1 t  yoe=  1 
1 1 id=0 
1  I iv  =  0 
1 1 ke=2 

1  I  n  1  =0 
1 ln2  =  0 
1 ln}  =  0 
1 ln<l=0 
1 1nS=0 
neoun=0 
1 ldof =6 

nicui 
qtvoe=*G' 

1 1 eonf =0 
1  I  c  i  o  =  0 
LSPC1=0 
LSPC2=0 
LSPC3=0 
LSPCasO 
LSPC5=0 
LSPCb-0 

INITIALIZATION  DONE 

k j  =  j  nobed- i nn i 

k  row  =  k  j -2 

do  200  jsl, k  row 
C  looo  on  new  nb.  of  row 

c 

nc  ounsncount lA(jnnitj-l) 

11otslA(jnnifj-l)-2 
do  100  1=1. Hot 
1 1  i d= 1 1 i d ♦  1 

xlloe=A<neoun*?tjxloct2*l) 
yl 1oes4(neoun*2»jxloet2* 1 ♦  1  ) 
e  skio  first  node  of  row 
e 

* l 1 oc-0 . 

if  (  4 ( j ve 1 o“2) • ne . 0 , )  then 
reouro=t./4(jve1o*2) 
thetas«lloe/rcouro 
xlloesreoure*sin<th«ta) 
tlloe=reourb*eos(th*ta)”rcourb 
enai  f 
c 

if  (a(jvelo”l).ne.0,)  then 
reouro=l ./a( j  ve  1  o •! ) 
theta =v 11 oe/reouro 
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vl  1ocsrcouro»sin(thefa) 
rl1ocsrcouro*eos(thet»)-reourb 
endi  f 


if((a(jve1o-l).ne.0.),and.(a(jve1o-2).ne.0.))  then 
nr i te( 1 1 u» 70 ) 

format (' error > tvo  curvatures  are  non  tero') 
endi  f 

ready  to  disolav  oacket 

*r»te(llu#30)  11tyoe*11id#11iv.l1kC/11nl,lln2,1lnJ,1ln4,11n5 
f ormat ( i 2«  8 i 8 ) 

-rited  lu.811  «  1  1  OC  <  v  1  1  OC  »  Z  1  1  OC 
format ( 3e  t  d . 0 ) 

vri te( 1 1u.32)  11 ief .dtyoe. 1 1  do f , 1 1  eon f , 1 1 c i d, LSPC 1 , LSPC2 , 
1  LSPC 3  LSPC4.LSPCS.LSPC6 

format(It#al#3f8,2a*6il) 

CSnt i "ue 
Conti nue 
ret  urn 

end 


subroutine  endneu 
Comnon/franie/nfntot#  1  1u 
1  1  tvoe=<»P 
1 1 i d*0 
1 1 ivru 
11kc*1 

vrite(11u>80)  11tyoe»11id»11iy,llke 
f  omat  ( i  2«  3i  8) 
ret  urn 
end 
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nnnnnoo  non  n  oon 


c 


SUBROUTINE  N£U0ISP(0£FL, VELO) 

01MENSI0N  4(1) 

COMMON  14(1) 

EOU I  VALENCE  (14(1),  4(1)  ) 

DIMENSION  OEFL ( 1 ) » V£L0( 1 ) 

COMMON  /CP4R4/  NSTEPS  ,  NSEG  ,  NENO  ,  NSHEET  ,  N28D  ,  N3BDI  , 

1  N3802,  INTRVL,  OELT,  NMTOT ,  NJOIN,  NRELAX,  4LPM4 

S  LEN  SB  X 

COMMON  /  SSIZE  /  IBG( 1 ) , N0J , NELTOT , N I  BO. NLOAO, NBRECT , 

1  N80U40,  ISHEET,  NPRPTS,  NSPTS,  NSTPPTS,  NVPTS,  NHPTS, 

2  NSSTrP,  NNJ,  NNTOT,  LGOSP,  LIQUD,  LBC4LC 

COMMON  /ST4B  /  I8T( 1) . JSSIZE, JSPAR, JVELO, JSTRE.JXMAS, JIELM, JBmaT, 

1  JLlBO. JtOOP, JPRET.JLHlS, JSTRN, JFORC, JXLOC, JNQI, JNNI, JNQ8EG, 

*  JLSIOE. JIELMCL, JSTIF, JOEFL, JFORLG, 

2  JIFPAR, JFLPAR.JXCOORO.JYCOORO, JOEL  TAX, JOEL TAY.JVMA.JSEFX, 

3  JFLUFR ,  JP9INC,  JVELR4 0 ,  JGENFR , JPRES , JCSEP 
COMMON/CIO/  NIN,N0UT,NTHIST,NC0RT,MC0RT,NTPL0T,NVMA,ITITlE(20) 
C0MM0N/FRANF/NFNT0T,LLU.LLN,LLS,NDPLTS,NSPLTS 

COMMON/ T I TRE/NTITL£(80) 

COMMON/COECT/ISTEP.TIME 


REPORT 
BL»NC  2 
BLANC  3 
BLANC  a 

CPARA  2 
CPAR4  3 
MS PAR A  3 
SSIZE  2 
SSIZE  3 
SSIZE  « 
REPORT  8 
STAB  2 
STAB  3 
STAB  a 
STAB  5 
STAB  b 


INTEGER  NSU8T1(90).NSU8T2(80) 

KJsJNQBEG-JNNI 

NFNTOTsNNTOT-IA (JNNI ) -I  4 ( JNOBEG-1 )-2«(FJ-2) 

LLU*  l P 
LLNS18 
LLMslZ 
0EFM4XS0. 

MAXNOOsNFNTOT 

NWtOTHaS 

NOMAXSO 

LLIO=0 

NCOUN*0 

FROw=F J-2 

00  10  J*l ,80 
NSUBT 1 ( J ) s  0 
10  NSUBT2(J)S0 

PUT  the  time  STEP  ON  FIRST  SUBTITLE: 

CLOSE (UNI T*LLM) 

0PEN(UNIT=LLM,ST4TUSs'NEW ) 

MR  I T£(CLM, * )  TI«e 
CLOSE  (UNIT  =LLM ) 

RE40(LLM,99)  (NSUBTl(J),Jst,80) 

CLOSE (UNIT *LLM) 

OPENIUNITsLLM.STATUSs'OLO' ) 

OPEN ( UN tT=LLN, FORM® 'UNF0RM4 TTEO' .STATUS *’ NEW' ) 

FROM  IS  THE  NEM  NUM0 .  OF  ROMS  ,KJ  IS  THE  OLO 

NCOUN  COUNTS  N004L  PTS  IN  EPSA,  LLID  COUNTS  FOR  ACTUAL  «OOEL 

NUMR  COUNTS  NOOAL  PTS  IN  EPSA 

FIRST  LOOP  TO  GET  MAX.  OEFORMaTION  AND  NODE  NUMBER 

00  200  Js 1 , FROM 
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NC0JN=NC0UNyIA(JNNI*J-1  ) 

LLPT=IA(JNNI»J)-2 

C 

C  LLPT  IS  NU«8.  OF  POINTS  (FDR  PATRAN)  IN  EACH  ROW 

C  i*E  TA*£  THE  RIGID  BOOT  MOTIDN  OUT  9Y  SUBTRACTING  THE  V  AND  ft 

C  DISPLACEMENTS  AT  END9LDCK  AT  EACH  NODAL  POINT 

C  NCDUN*2  IS  THE  POINT  N8 .  (EPSA)  FDR  ENQ9LDCK 

C 

R9Y=DEFL(3*(NCDUN*2)-1) 

R8Z=0EFL( J*(NC0UN*2)) 

IF  (LIOUO.EQ.O)  THEN 
R8Y=0. 

R8Z=0. 

ENOIF 

00  100  L  =  1 » LLPT 

LLlO=LLIO*t 

NUMR=NCDUN»L*1 

OX=OEFL(3»NUMR-2) 

DY=DEFL(3*NUMR-1 )-RBY 
OZsOEFL(3*NU'»R)-RBZ 
C 

C  IS  SHELL  CURVED, CHANGE  COORDINATES 

C 

IF( (A(JVELO-l) .NE.O.) .OR. (A(JVELD-2) .NE.O. ) )  THEN 
CALL  CHCDDRD(NUMR,Dx,DY,DZ) 

ENOIF 

OD=AMAXl (A9S(0X) , ABS(DY) , A8S(DZ) ) 

IF(A8S(0EFMAX) .LT.DO)  THEN 

NOHAXSLLIO 

OEFMAXsOO 

IFCABSCAMINl (OX,DY,OZ)).EO.DO)  THEN 
OEF*AXs-DEF«AX 

ENOIF 

C 

C  he  CMECKEO  IF  DEFHAX  HAS  NEG. 

C 

ENOIF 

100  CONTINUE 

200  CONTINUE 

C  OK  FDR  FIRST  CARO 

WRITE(LLN)  (NTITLE(I) , I=! , 80 ) , NFN TOT , MAXNDO , DEFM AX , ND«A X , 
t  NHIOTH 

wrt  te(LLM»90 )  TITLE,NFNT0T,MAXNDD,DEFMAX,NDMAX,NMIDTH 

90  FDRHATCA, 218. £16.9,218) 

95  FORHAT(20At,2I8,Et6.9,2I8) 

HRITE(LLN)  (NSUBTl(I), 1=1,80) 

HRITECLLN)  (NSU8T2 ( I ) , I = l , 80 ) 

HRITE(LLH,9l ) 

91  FORMATCPINE’ ) 

99  forma t(80ai  ) 

LL 1 0  =  0 
NCDUN=0 
C 

C  SECOND  LOOP  TO  PICK  UP  DEFLECTION  AT  each  NODE  (OF  PATRAN  MOOED 

c 

00  900  Jsl.KRON 

NCDUN=NCDUN*IA(JNNI»J-t  ) 

LLPT=IA( JNNl»J)-2 
C 

C  take  OUT  RIGID  BODY  MOTION,  RBX.RBZ  DISPLACEMENTS  AT  END  BLOCKS 

C  ASIJMED  TO  REPRESENT  RIGID  BOOT  MOTIONS 


100 


c 


rby=oefl(3* cncoun*2)  -n 

S8ZsOEFL( 3*<NC0UN*2) ) 

IF  TLIQUO.EO.O)  THEN 
RBYsO. 

RBZxO. 

ENOIF 

C  IF  NO  FLUlO  OPTION  00  NOT  SUBTACT  RIGIO  300Y  CONTRIBUTION 
DO  500  L=1.LLPT 
LLI0=LLI0»1 
NUMRsNCOUNyL* 1 
0X=0EFL(3»NUMR-2) 

OY=OEFL(  5*NU«R-n-RBY 

OZ=OEFL(5»NU^R)-R8Z 

OZLOCsOZ 

VELZ=VEL0(3»NUMR) 

I F  < (A{ JvELO-1) .NE.O.  ) .OR. CA< JVELO-2) .NE.O.))  THEN 
CALL  CHCOORO(NUMR,DX,OY,DZ) 

ENOIF 

rtRITE(LLN)  LLIO.OX.OY.OZ.OZLOC,  VELZ 
*RITE(LLM,92)  LLI 0, OX , DY , OZ , OZLOC , VELZ 
300  CONTINUE 

«00  CONTINUE 

C  CLOSE  FILE  OPENEO  ON  UNIT  LLN 
CALL  CLOSE(LLN) 

92  FORMAT! I8.5E16. 9) 

RETURN 

ENO 

C 

c 

c 

SUBROUTINE  CHCOORO ( N , X , Y , Z ) 

0 1 MENS I  ON  A ( | ) 

COMMON  I A ( t ) 

EQUIVALENCE  ( I A  ( 1 ) ,  A ( | )  ) 

COMMON  /CPARA/  NSTEPS  ,  NBEG  ,  NENO  ,  NSHEET  ,  N2B0  ,  N3B01  , 

1  N3B02.  INTRVL,  BELT,  NHTOT,  NJOIN,  NRELAX ,  ALPHA 

S  LEN  SBX 

COMMON  /  SSIZE  /  IBG< 1 ) ,NQJ,NELT0T,N1B0,NL0A0,NBRECT, 

1  NBQUAO,  ISHEET,  NPRPTS,  NSPTS,  NSTRPTS,  NVPTS,  NHPTS, 

2  NSSTYP,  NNJ,  NNTOT ,  LGOSP,  LIO'JO,  LBC ALC 
C 

COMMON  /STAB  /  IBT( 1) , JSSIZE, JSPAR, JVELO,JSTRE,JX«AS,JIELM,JBMATt 
l  JL1B0. JLOOP,JPRET,JLHIS,JSTRN, JFORC, JXLOC,JNQI, JNNI,JNQBEG, 

1  JLSIOE,JIELMCL,JSTIF, JOEFL, JFORLG, 

2  JIFPAR,JFLPAR, jxcoopo.jycooro, joeltax, JOELTAY, JVMA.JSEFX, 

3  JFLUFR,  JPRINC,  JVELRAO,  JGENF R , JPRE S , JCSEP 

common/cio/  nin,nout,nthist , ncort ,mcort » nt plot , nvma, i title (20) 

COMMON/FRANK/NFNTOT,LLU,LLN,LLS,NOPLTSfNSPLTS 

DIMENSION  XMAT (3,3) 

00  10  1=1,3 
00  10  J* 1 , 3 
XMAT(I,J)x0. 

10  CONTINUE 
C 

C  THIS  ROUTINE  CHANGES  THE  OISPL.  OF  NOOAL  PT.  N ( FOR  EPSA) 

C  IN  A  GLOBAL  RECTANGULAR  SYSTEM 

C 

C  IF  CURVATURE  IN  Y  OIR.  IS  NON  ZERO  , CALCULATE 

C  The  ROTATION  matrix  AT  each  POINT 

C 


REPORT  2 
BLANC  2 
BLANC  3 
BLANC  « 
CPARA  2 
CPARA  3 
MSPARA  3 
SSIZE  2 
SSIZE  3 
SSIZE  « 
REPORT  8 
STAB  2 
STAB  3 
STAB  9 
STAB  5 
STAB  b 


IF  (At  'y£L0-1  )  .'(£  .  I).  1  TH£'t 

aca’javsi ./»  (  JVELJ-i ) 

TH£YAI*A(Jxl0C*2*N-1  1/aCOOS* 

xmaui,  i  j*i. 
xmat(2,2)*C0S(ThETAI  ) 

XMAT (2, J)sSIN(Th£TA1 ) 

XMAT(3,3)*<MAT(2,2) 

XKAT(J,21S-1..XMAT(2. si 
CALL  P«00( X«AT , X , r , l ) 

END  If 
C 

c  IF  CURVATURE  IN  *  DIRECTION  IS  NON  ZERO: 

C 

IF  (A(jv£L0-2).N£.0.)  THEN 
DO  20  1*1.3 
00  20  J*1  .  3 
XMAT(I,J)*0. 

20  CONTINUE 

RCOURX* 1 ./A( JVELO-2) 

THETA2=A(JXL0C*2,N“2l /RCOURx 
XMAT(2.2)*1. 

XMAT ( 1 , 1 ) *COS ( Th£TA2 l 
XMAT(3,3)*XMAT(l,l) 

XMAT (3, l  )=-l ,*SIN(TH£IA2) 

XMAT(1,3)=SIN(ThETA2) 

CALL  PROOCXMAT.X.Y.Z) 

ENOIF 
RETURN 
END 
C 

c 

SUBROUT  I N£. PROD (XWAT.X.Y.Z) 

DIMENSION  XMAT (3,3) 

C 

C  THIS  ROUTINE  DOES  The  MATRIX  PROOUCT  XMAT.(X,Y,Z) 

c 

Xl=X 
Y  1  *  Y 
21  =  2 

XsXMAT (1 , 1 )  *  X  l tXMA  T ( 1 ,2) *Yl*XMAT ( l , 3) *21 
YaXMAT(2, 1}»X1 fXMAT(2,2)»YI t  X M A T  ( 2 , 3  )  *Z  l 

Z*xmaT(3,1)»xi»xmat(3,2)»Y1yxmaT(3,J)*Z1 
return 

ENO 
C 

c 

SUBROUTINE  PLOTOISP (OEFL, VELO ) 

OIMENSION  OEFL ( I ) » VELO ( l ) 

COMMON/COELT/ISTEP.TIME 

COMMON  /CPARA/  NSTEPS  ,  NBEG  ,  NENO  ,  NSHEET  ,  N2BO  ,  N3BOI  , 

1  N3B02,  INTRVL,  OELT ,  NHTOT,  NJOIN,  NRELAX,  ALPHA 

S  LEN  SBX 

COMMON/FRANK/NFNTOT,LLU,LLN,LLS,NOPLTS,NSPLTS 

c 

C  CHECKS  IF  OUTPUT  FOR  PATRAN  IS  REQUESTED  BY  USER  IF  YES  CALL 

C  NEUOISP 

C 

IF  (TIME, EQ. OELT)  XOISP*l 

timei*float(nsteps)/float(noplts) 

TIM£C*K0ISP*TIM£I 

ITIM£aJNlNT(TIMEC) 


CPARA  2 
CPARA  3 
MSPARA  3 
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TiMEcsT[Me/(in«e*oEi.n 

EVOTIMsNSTEPS*OELT 

TESTS A9S( TIMEC-JNINTCTIMECI ) 

IF (<  TEST. LE. 0.0000 1  ) .OP. <Tr*E.EQ.ENOTIM) )  THEN 
call  neudisp(OEfl,v£lO) 

KOISPsKOISP+1 

ENOIF 

RETURN 

ENO 


SU8R0UTINE  PLOTSTRE(STRE) 

OIMENSION  STRE(l) 

C0MM0N/C0ELT/IST£P,TIME 

COMMON  /CP AR A /  NSTEPS  ,  N8EG  ,  NEND  ,  NSHEET  ,  N2B0  ,  N3601  ,  CPARA  2 

1  N3BD2,  INTRVL,  OELT,  NHTOT,  N JO IN,  NRELAX,  ALPHA  CPARA  3 

S  LEN  SB*  MSPARA  J 

COMMON/FRANK/NFNTOT,LLU,LLN,LLS,NOPLTS,NSPLTS 

C  CHECKS  IF  OUTPUT  FOR  PATRAN  IS  REOUESTEO  ST  USER  IF  YES  CALL 

C  NEUOISP 

c 

IF  (TIME. EQ. OELT)  KSTPEsl 

TIME  I SFLOATI NSTEPS )/FLOAT(NSPLTS) 

TIMEC=KSTRE«TIMEI 

ITIME=JNINT(TIMEC) 

TIMEC=TIME/(ITI«E*OELT) 

ENOTIM=NSTEPS*OELT 

TESTsABS(TlMEC-JNINT(TIMEC)) 

IF( (TEST. LE.O. 0000 1). OR. (TI«£.EQ.£NOTIM))  THEN 
CALL  N£ JSTRE ( STRE ) 

KSTRE*KSTRE*1 

ENOIF 

RETURN 

ENO 

C 

C 


C 


SUBROUTINE  NEUSTRE(STRE) 

OIMENSION  A (  1 ) 

COMMON  I A ( 1 ) 

EQUIVALENCE  (Ia( 1 )  ,  A(l)  ) 

OIMENSION  OISP( 10) ,STRA(10) ,OSTRE< I0),0S(5) , VEL(2) 

OIMENSION  STRE ( 1 ) 

COMMON  /CPARA/  NSTEPS  ,  N9EG  ,  NENO  ,  NSHEET  ,  N2B0  ,  N3B01  , 

I  NJB02.  INTRVL,  OELT,  NHTOT ,  NJOIN,  NRELAX,  ALPHA 
S  LEN  SBX 

COMMON  /  SSIZE  /  I8G< 1),NQJ,NELT0T,N1B0,NLDA0,NBRECT, 

1  NBOUAO,  ISHEET,  NPRPTS,  NSPTS,  NSTRPTS,  NVPTS,  NHPTS, 

2  nssttp,  nnj,  nntot,  lgosp,  liquo,  lbcalc 

COMMON  /STAB  /  IBT(1),JSSI2E, JSP4R,JVEL0, JSTRE, JXMAS,JIELM,J8MAT, 

1  JL180, JLOOP, JPRET, JLHIS, JSTRN, JFORC.JXLOC, JNQI, JNNI,JNQBEG, 

$  JLSIOE, JIELMCL, JSTIF, JOEFL, JFORLG, 

2  JIFPAR, JFLPAR, JXCOORO, JYCOORU, JOEL  TAX, JOEL  TAT, JVMA, JSEFX, 

J  JFLUFR ,  JPRINC,  JVELRAO,  JGENFR , JPRES, JCSEP 

common/cio/  nin,nout,nthist,ncort,mcort,ntplot,nvma, ititlE(20) 

COMMON/FRANK/NFNTOT,LLU,LLN,LLS,NOPLTS,NSPLTS 

C0MM0N/TITRE/NTITLE(80) 

COMMON/COELT/ISTEP.TIME 


REPORT  2 

blanc  2 

BLANC  3 
BLANC  a 


CPARA  2 
CPARA  3 
MSPARA  3 
SSIZE  2 
SSIZE  3 
SSIZE  a 

REPORT  8 
STAB  2 
STAB  3 

stab  a 

STAB  5 
STAB  t> 
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o  o  o  o  o  n  o 


COMMON/SPARA/HEGO  )  »£,GNU  RHO.SGYLO,  THICK  ,  CUR  V  (  2  ) 

INTEGER  NSU9T1  (80)  ,.NSU9T2(80  ) 

NwIOTHsll 
LLS- 1 6 
LLZ=1 5 

00  10  J=l,80 
10  NSU8T2(J)=0 

C  TAKg  CARE  OF  The  2N0  TITLE  ,  INOICAT.  TIME  step 

c 

CLOSE (UNI TsLLZ) 

OPEN  (UN  I  TsLLZ.  STATUS3 'NEW  ) 

WRITE (LLZ  .  * )  TIME 
CLOSE(UNITsLLZ) 

REAO(LLZ.SOl)  (NSU8Tim,I  =  l,80) 

CLOSE(UNlTsLLZ) 

OPEN(UNIT=LLZ»STATUS='OLO* ) 

801  FORHAT(SOAl) 

OPEN (UNI T=LLS, FORHs 'UNFORMATTED1 , STATUS: • NEW ) 

C 

C 

C  WE  ARE  GOING  TO  OISPLAY  THE  FIRST  THREE  RECORDS  (TITLES) 

C 

WRITE (LLS)  (NTITLE(I) ,l3t ,80),NWIDTH 
WRITE(LLS)  (NSU9T1(I),I=1,80) 

WRITE(LLS)  (NSU8T2(I),lsl,80) 

WE  are  GOING  TO  SCAN  ALL  ELEMENTS  AS  IN  SHCONN, 

GET  THE  STRESSES  PERFORM  ROTATION  IF  NECESSARY 
LLELl  NU«.  OF  rows  ELTS:  LLIOlELT  number 
NPREV:  NOOE  NUM, (FOR  EPSA)  OF  LAST  POINT  OF  ROW  JUST 
8EL0W  THE  ELEMENT  ROW  K;  LLROw:  NUM.  OF  ELTS  IN  EACH  ROW 

NPREVSO 

C  IA  ELTS  ♦! (FOR  NOOES)  *2  ARTIFICIAL  NOOES  FOR  NPREV 

NSHAPEua 

llel=jnni-jnqi 
00  200  Kst.LLEL 
C 

LLR0W=IA(JNQI»K-1) 

NPR£VsNPREV*LLR0W»3 

C 

00  100  J=l, LLROW 
LLl0*lA(JN08£G»K-l)»J-l 
03(1 )sSTRE(9»(LLI0-t ) M ) 

OS(2)sSTREO»(LLlO-l)»2) 

OS (3  )30 . 

C  OS ( u . 5 )  WILL  SE  THE  STRESSES  IN  LOCAL  COORO. 

OS(a)=OS( 1 ) 

OS  ( 5 ) *03 ( 2 ) 

C  VON  MI3ES  AT  TOP  ANO  BOTTOM  OF  SHELL: 

T0XYSSTRE(R»(LL10-1 )♦!) 

OSXsSTREO*(LLIO-l)t«) 

03YaSTRE(R*(LLlO-l)fS) 

OSXsOSX*6. /THICK 
03Ys0SY»6. /THICK 
03Xls0S(l)»0SX 
03Yt*0S(2) »OSY 
0SX2*0S( 1 J-OSX 
03Y2*0S(2)-0SY 
C 

C  VON  MISES  CALCULATIONS  MODIFIED  AFTERWAROS  FOR  MEMB. *8EN0ING  STRESS 


o  r»  n  o 


NO  SMEAR  at  top  of  smell, get  stresses  at  TOP  ANO  BOTTOM  ANO  take 
max,  value 

top: 

SIGM1=0SX1 
SIGM2SOSY1 

VONMlsSORTCSIGMl«SIGMl-SIGMl «S IGM2 »S IGM2.S IGM2 ) 

C  90TT0M; 

SIGMHOSX2 
S  IGM2=OSY  2 

VONM2=SQRTCSIGMl •SIGM1-SIGM1 «S IGM2  *S IGM2»S IGM2 ) 

C 

VONM=AMAXl CV0NM1,V0NM2) 

C  O.K.  FOR  VON  RISES  CVONM) 

C  IS  THE  SHELL  CURVEO,  NLELI’S  ARE  NOOES  SURROUNOING  ELT.  (EPSA) 

c 

IF((A(  JVELO-1)  .NE.O.)  .3R.CACJVELQ-2)  .NE.O.n  THEN 
NL£L1=JaN?R£V*! 

NLEL2=J*NPREV*2 
NLEl3=J*1*NPR£V*LLROwm *2 
NLEL«*J»NPREV*LLR0Na1*2 

XsCA(JXLOC*2*NLEL1-2)  +A  t  JXL0C«2»NLEL2 -2  ))/<». 
X=x*(A(JXLOC»2*  NLEL3-2 ) »A  t  JXLOC ♦2*nlELV-2 ) )  /« . 

Y  =  CACJXL0C»2*NL£LI-I  )  *AC  JXL0C*2«NLEL2-1  n/a. 
YsY*(A(JXLOC*2*NLEL3-t) ♦ A  C  JXLOC  *2*NLELa- I ) )/R. 

CALL  CHELEMCX.Y.OS] 

ENOIF 

C 

C  REAOY  TO  OISPLAY  RECORO 

C 

hRITECLLS)  LLIO.NSHAPE, COISPCl) ,1=1, 10) , 

1  (STRA(I) ,1=1, 10! .0STREC1 ) ,0SC5) .03(6), 

2  OSTREC«),COS(I).I=t.3)» 

3  COSTREC II , 1=8, 10) , VONM 

C  hRITE(15,30)  LLI0.N3HAPE, (03(11,1=1,5), VONM 

100  CONTINUE 

200  CONTINUE 

CALL  CLOSECLLS) 

80  FORMAT (2I0,6E15.8) 

RETURN 

ENO 


C 


SUBROUTINE  CHELEMCX.Y.OS) 

DIMENSION  a(1) 

COMMON  I A ( 1 ) 

EQUIVALENCE  CIAC1),  A(l)  ) 

COMMON  /CPARA/  N3TEPS  .  N3EG  »  NENO  ,  NSHEET  ,  N2B0  ,  N3B01  , 

1  N3B02,  INTRVL.  OELT ,  NHTOT ,  NJOIN,  NRELAX ,  ALPHA 

S  LEN  SBX 

COMMON  /  SSIZE  /  I8G(1),NQJ,NELT0T,N1B0,NL0A0,N6RECT, 

1  NBGUAO ,  IJHEET,  NPRPTS,  N3PT3 ,  NSTRPTS,  NVPTS,  NHPTS, 

2  NSSTYP.  NNJ,  NNTOT,  LGOSP,  UOUQ,  LBCALC 

COMMON  /STAB  /  IBT ( J) , J5SIZE, JSP AR , JVELO, J3TRE, JXMAS, JIELM, JBMAT , 

1  JL1BQ,JL00P,JPRET, JLHIS, JSTRN, JFORC, JXLQC,JN0I, JNNI,JNOBEG, 

J  JLSIQE ,  J IELMCL ,JSTIF,JOEFL,JFORLG, 

2  JIFPAR, JFLPAR,JXCOORO,JYCOORO, JOEL  TAX, JOELTA Y , JVMA , JsEFX , 

3  JFLUFR,  JPRINC,  JVELRAO ,  JGENFR , JPRES, JCSEP 
DIMENSION  XMAT (3, 3) 

OIMENSIQN  OS ( 3) 
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o  o  o  o  o  o  o 


10 


oo  to  in,  5 
00  10  J  =  1.5 
xurti.jjso, 

CONTINUE 

This  SOUTINE  CHANGES  TriE  STRESSE  IN  AN  ELEMENT  INTO  A 
IN  A  GLOBAL  RECTANGULAR  SYSTEM  .GIVEN  LOCAL  COORO.  OF 
CENTRO  10  x.Y 


IF  CURVATURE  IN  Y  OIR.  IS  NON  ZERO  .CALCULATE 
THE  ROTATION  MATRIX  AT  EACH  POINT 
IF  (A(JvELO-t).NE.O.)  THEN 
RC0I)RY  =  1  ./A(JVEL0-1  ) 

THETA1=Y/RC0URY 
XMAT(l, |)sl  . 

XMAT(2,2)sC0S(THETAl) 

XMAT(2,J)sSIN(THETA1) 

XMAT(3,3)=XMAT(2,2) 

XMAT( J,2)s-t.»XMAT(2, J) 

CALL  PROO(XMAT,DS( 1 ),0S(2),0S(3)) 

ENOIF 
C 

C  IF  CURVATURE  IN  X  OIRECTION  IS  NON  ZERO: 

IF  ( A ( JVELO-2)  .NE.O.)  THEN 
00  20  1=1,3 
00  20  J  =  l,3 
XMAT(I,J)=0. 

20  CONTINUE 

RCOURXsi ,/A( JVELO-2) 

THETA2=X/RC0URX 

xmat (2 , 2 ) s 1  . 

XMAT ( 1 , 1 ) =C0 S I THET  A2) 

XMAT(J,3)sXMAT(l,l) 

XMAT (3, 1 )S-1 .*SIN(THETA2) 

XMATU  ,3)=SIN(ThETA2) 

CALL  PROO(XMAT,OS(! ) ,0S f 2 ) , OS ( 3 ) ) 

ENOIF 

RETURN 

ENO 
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